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CONTROL OF CHAOS IN LASERS 
BY FEEDBACK AND NONFEEDBACK METHODS 


PIERRE GLORIEUX 
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Centre d’Etudes et Recherches sur les Lasers et Applications, 
Universite des Sciences et Technologies de Lille, 
F-59655 Villeneuve d’Ascq Cedex, France 

Received July 31, 1997; Revised November 25, 1997 


Analytical approaches of the theory of chaos control in class B lasers are presented in the 
cases of continuous delayed feedback and of subharmonic modulation. They are compared with 
experimental findings on a CO 2 laser and a Nd-doped fiber laser with modulated loss and 
pump respectively. In both cases, analytical theory allows one to predict the shift of the first 
period-doubling bifurcations. Numerical simulations show that subharmonic modulation may 
induce shifts, crises and new attractors in a laser with modulated parameters and that its phase 
relative to the main modulation allows one to control the laser dynamics. These results agree 
well with the experimental observations on control of chaos in CO 2 and fiber lasers. 


1. Introduction 

It is now conventional to separate among the 
methods used to control chaos those which rely on 
feedback from those based on nonfeedback tech¬ 
niques. In the former some information from the 
chaotic system is extracted and used to design a 
weak correction signal applied to some suitable con¬ 
trol parameter. Following the seminal work of Ott, 
Grebogi and Yorke [Ott et al, 1990], several al¬ 
gorithms have been proposed to control chaos by 
stabilizing a system on one of its unstable peri¬ 
odic orbits and applied to various systems includ¬ 
ing lasers [Roy et al, 1994: Bielawski et al, 1993]. 
In nonfeedback control, a small periodic perturba¬ 
tion or parametric modulation is shown to drasti¬ 
cally alter the dynamics of the system under con¬ 
sideration [Bryant & Wiesenfeld, 1986] and this was 
applied to both monomode [Meucci et al, 1994; 
Chizhevsky et al, 1986; Gavrielides et al, 1985] 
and multimode lasers [Otsuka et al, 1997]. These 
two approaches are of completely different natures, 
but they share the concept of controlling dynam¬ 
ics through the use of small perturbations. The 


weakness of these perturbations suggests the exis¬ 
tence of small parameters in both situations and 
therefore there is the possibility of applying pertur¬ 
bation methods to obtain analytical properties of 
the systems subjected to these control techniques 
and possibly to compare them with experimental 
findings. 

In this paper we report on advances in the 
dynamics of class B lasers in presence of (i) de¬ 
layed feedback [Bielawski et al, 1994; Erneux 
et al., 1995] and (ii) subharmonic modulation 
[Newell et al, 1997; Dangoisse et al, 1997]. These 
works show that singular perturbation methods and 
more precisely multiple timescale analysis is able to 
predict how the dynamics are altered by the control 
techniques and in particular to predict if bifurcation 
points are shifted by the use of the control tech¬ 
niques mentioned above and to describe how this is 
done. All these predictions are compared with ex¬ 
perimental results on chaos control of class B lasers. 
This point of view is different from that developed 
in [Chizhevsky et al., 1986] where the emphasis is 
on the fact that the nonfeedback control of chaos 
creates new attractors, a point which was already 
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discussed for a different system by Gavrielides et al. 
[1985]. 

Control of chaos by delayed feedback was 
proposed by Pyragas [Pyragas, 1992; Pyragas & 
Tamasevicius, 1993] and applied in the laser do¬ 
main to loss modulated CO 2 lasers [Bielawski 
et al, 1994; Erneux et al, 1995]. Chaos may 
also be controlled by adding a small modulation 
which is resonant or nearly resonant with a subhar¬ 
monic of the fundamental system frequency [Meucci 
et al, 1994; Chizhevsky et al, 1986; Braiman & 
Goldhirsch, 1991; Colet & Braiman, 1996; Ciofini 
et al, 1995; Vohra et al, 1995]. This technique 
was first applied to a magnetostrictive ribbon and 
soon after to lasers including CO 2 [Meucci et al, 
1994; Chizhevsky et al, 1986], microchip [Otsuka 
et al, 1997] and fiber [Dangoisse et al, 1997] lasers. 
All these belong to the so-called class B lasers in 
which the relaxation time of the electric polariza¬ 
tion is much shorter than the photon cavity life¬ 
time which is itself much smaller than the popula¬ 
tion inversion relaxation time. Class B lasers have 
been used as test systems in this approach because 
they display well characterized chaos when they are 
subjected to periodic modulation of one of their pa¬ 
rameters [Ivanov et al, 1982; Arecchi et al, 1982; 
Tredicce et al, 1986; Chen et al, 1985]. Moreover 
single-mode class B lasers are efficiently described 
by two-variable models suitable for fully analytical 
treatments. More specifically, doped fiber and CO 2 
lasers have been used in the corresponding experi¬ 
ments. These particular lasers have been chosen for 
technical reasons linked to timescales and accessi¬ 
bility of parameters required for control. In fact the 
fiber laser is not a single mode laser but the stan¬ 
dard class B model correctly describes its dynamics 
within some parameter range discussed in previous 
papers [Derozier et al, 1992]. 

Chaos in class B lasers has long been observed 
when a suitable parameter is modulated at a fre¬ 
quency in the domain of the free running relax¬ 
ation oscillations. In CO 2 lasers, chaos has been 
observed when the pump parameter [Biswas et al, 
1987], the cavity loss [Arecchi et al, 1982; Tredicce 
et al, 1986] or the cavity detuning [Midavaine 
et al, 1985] have been sinusoidally modulated. The 
relative efficiency of these modulations has been 
discussed [Khanin, 1995] but technical arguments 
eventually determine which is the most efficient for 
a specific device. In doped fiber and YAG lasers, 
pump modulation has almost exclusively been used 
for the same technical reasons. 


A single mode class B laser is essentially a 
two-dimensional dynamical system which is well de¬ 
scribed by the model 

f t =2I[AD-l-k{t)] 

§ -Mi-JXi + fl] 

The variables I and D are the laser intensity and 
population inversion respectively. Time t is mea¬ 
sured in units of the cavity lifetime. A is the pump 
parameter and 7 is the ratio of the population in¬ 
version rate to the cavity damping rate. k(t ) is the 
additional loss term containing the modulation and 
possibly feedback terms. The parameter 7 « 10 -3 
is small for any common class B laser which in¬ 
cludes also semiconductor in addition to CO 2 and 
solid state lasers. By removing this factor from 
the right hand side of the second equation, one ob¬ 
tains equations for the deviations x and y from the 
zero-intensity steady-state (I, D) = (Jo, Do) where 
Iq — A — 1 and D 0 = 1/A and introducing new vari¬ 
ables and parameters as defined in [Erneux et al, 
1995] eventually leads to the basic equations for the 
free-running class B laser 

f = -y — ex[l + (A - 1)(1 + y)\ 
ts ={ 1+y)X 

where s is the new time scale s = (27 0 7) 1 ^ 2 f and 
e=y/2 7 (A-l). 

In presence of loss modulation with an adimen- 
sional efficiency M and a delayed feedback with ef¬ 
ficiency f3, the second equation should be changed 
into 

dv 

j- ? = (l + y)[x - P(y - y(s - v)) -M cos as] 

where v is the time delay of the feedback loop and 
<7 the reduced modulation frequency. 

If the pump is modulated, e.g. by a two-tone 
signal as in the second half of this paper where 

A = Ao(l + r\ cos u\t + r 2 cos+ ip)). 

The first equation should be altered to account for 
this modulation and reads then as 

fjqf 

Ts =- y -ex[l + (Ao-l)(l + y)} 

+ Si COS tJ\S + S 2 COs(<72S + p) 
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where Si and cr* (i = 1, 2) are the reduced modu¬ 
lation amplitude and frequency respectively for the 
two modulations. The phase </? is relevant only when 
<7 i/<72 is a rational number. 

Suitable methods are developed in each case 
with the aim of determining the position of the 
first bifurcation, i.e. the value of the parameters for 
which the period 1 solution destabilizes. This gives 
the parameter range in which the control is efficient 
in stabilizing the period 1 orbit. 


2. Analytical Analysis 

of the Delayed Feedback 

For the delayed feedback, asymptotic analysis of the 
equations has been carried out using e as a small 
parameter. The Oth order (e = j9 = 0 , M < 1 ) 
solution of the equations of the linear problem is a 
combination of 2 - 7 r-periodic and 2-n /cr-periodic func¬ 
tions of the form 

M 

X = 2.RcOs(s + 9) + -2 C0S<7S 

y = 2R sin(s + 9) + \ s j n as 

where R and 9 are unknown amplitude and phase 
respectively. We wish to determine the bifurcation 
diagram of the 2it/a and Air/a periodic solutions. 
Therefore we concentrate on a « 2 and consider the 
above equation as our leading approximation. The 
bifurcation analysis consists of studying the ampli¬ 
tude R = 0(e 1 / 2 ) and phase 9 = 0(1) as functions 
of the slow time es. They satisfy 

- 7 ^ = — ( ~eA + ft) R + \mRcos 29 
as V 2 / 6 

R~ = 1(2 - ff )R - -MRsm29 - -R 3 . 
ds 2 v ' 6 6 

Steady-state solutions correspond to periodic solu¬ 
tions of the original equations. Period 1 (R = 0, 
9 arbitrary) and period 2 solutions are determined 
and their linear stability analysis is performed. Let 
M = Ml = 13 (2 — a) | be the limit point from which 
two distinct branches of steady-states 9 = 9±(M) 
emerge. If M < Ml, there is no steady state for 9, 
which means that the laser is unlocked in phase. 

The stability condition of the period 1 solution 
is 

M l < M < M* 


with 


M* = 


Ml + 36[±eA + p 2 


1/2 


= Z[(eA + 2(3) 2 + (2 — cr) 2 ] 1/2 . 


If /3 > f3 c = — 1/2eA, a stable period-doubling 
transition is observed at M — M*(J3) where M* is 
defined above. If M < M *, the intensity of the laser 
is 


1-1 0 


JqM sin2s 


where M = 0{e). If M > M*, the leading approx¬ 
imation for the intensity is 


I — Io ~ — 2/o-Rsin(s + 9) 

where R = 0({M - M*) 1 / 2 ). 

If (3 < f3 c , a period-doubling transition is still 
possible but all periodic solutions are unstable. 
Quasiperiodic oscillations are then expected. 

The period 2 solution which corresponds to the 
R / 0 solution is stable when 

R 2 > 3(2 — a) if a <2. 

This period 2 solution is subcritical (M < M*) if 
a < 2 which implies instability for small R. Sim¬ 
ilarly the period 2 solution is supercritical (M > 
M*) if cr > 2 and is always stable because the above 
relation is always verified. 


3. Numerical and Experimental 

Investigation of Delayed Feedback 

The corresponding experiments were carried out on 
a CO 2 laser in presence of sinusoidal loss modula¬ 
tion. An electro-optic modulator inserted inside the 
laser cavity allows modulation of the cavity losses. 
The feedback may be easily applied by just adding 
the small correction signal to the main driving volt¬ 
age applied to the modulator. The delay is realized 
thanks to an optical fiber delay line which is 600 
meters long. The CO 2 laser intensity is detected by 
a photovoltaic detector and fed into a laser diode 
power supply. Part of the diode laser emission is 
directly detected and another part is sent into the 
fiber, allowing one to build a correction signal pro¬ 
portional to I(t) — I(t — t). The balance between 
the two channels is obtained by inserting an atten¬ 
uator in front of one detector and adjusting it so 
that the difference 1(f) — I(t - r) in presence of a 
T(= t ) periodically modulated signal is null. 
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Figure 1 shows the results of a typical experi¬ 
ment. The upper part is a bifurcation diagram of 
the modulated laser in the absence of feedback and 
is given as a reference. The lower part is the di¬ 
agram obtained with the same laser parameters in 
presence of delayed feedback. This feedback stabi¬ 
lizes the period 1 solution by shifting the period¬ 
doubling bifurcation point. Measurements of the 
shift magnitude versus the feedback coefficient are 
reported on Fig. 2. The experiment compares well 
with the numerical and analytical predictions in 
two points: (i) the bifurcation point shifts contin¬ 
uously to larger modulation values as the feedback 



t -—r- i — 1 ——ir^ 

5 10 15 20 

modulation amplitude (V) 

Fig. 1. Comparison of bifurcation diagrams (a) in the ab¬ 
sence and (b) in the presence of delayed feedback for a CO 2 
laser with modulated loss. 
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amplitude of the feedback loop (arb. units) 

Fig. 2. Shift of the first period-doubling bifurcation versus 
amplitude of delayed feedback as observed on a modulated 
CO 2 laser with delayed feedback. 


is increased from zero, and (ii) no period-doubling 
bifurcation is observed for negative values of the 
feedback. The analysis predicted a possible range 
of stabilization with negative feedback but it is too 
small to be observed here. The flat part of this curve 
in the region of small feedback is most probably due 
to discrepancies between the expected and the ob¬ 
served delays, which implies a loss of efficiency for 
low amplitude feedback. 


4. Control of Chaos by the Phase 
of Subharmonic Modulation. 
Analytic Results 


For this system, Erneux et al. obtained a mapping 
by using a method based on matched asymptotic ex¬ 
pansions and developed in [Newell et al ., 1997]. The 
pump modulated laser mapping relates the ampli¬ 
tude x n and the time interval s n between successive 
spikes. In the limit e -» 0, the period s n+ i — s n and 
the change in amplitude x n+ \ — x n are linked by 


25\ . . . 

s n+ i -s n = - 2 x n -I-sm(wis n ) 

LO\ 


+ — sin(w2Sn + V?) 

Ul 2 

x n+ i -x n = —-[sin(u;is n +i) + sin(wis„)] 

Wi 


5 2 

-[sin(w 2 s n + <p) + sin(w 2 Sn + V 3 )] 

0)2 

where the parameters describe the modulation of 
the pump parameter which writes in reduced units 
Si sin(wis n ) + 82 sin(w 2 s„ + (p ). 
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A fixed point analysis shows that if 82 = 0, 
the period 1 solution undergoes a period-doubling 
bifurcation at = 1. Erneux et al. investigated 
the bifurcation diagram of the periodic states in 
the vicinity of this point and for small values of 82 - 
They showed that if ip / 7r/2, the period 2 pitchfork 
bifurcation is destroyed and replaced by a smooth 
transition branch bounded by a fixed point. This 
limit point is located at 


tfi = 1 + 


6<5 2 sin 



In other words the standard bifurcation is then 
replaced by an imperfect bifurcation with two 
branches of solutions. For <p = tt/ 2, the imper¬ 
fection term disappears and the two-tone modu¬ 
lation problem reduces to the perfect bifurcation 
case. They also showed that when 82 is small, 
the period-doubling bifurcation is located approx¬ 
imately in 1 + 8 $. 


5. Experimental and Numerical 

Investigations of Control by the 

Phase of Subharmonic Modulation 

The influence of the phase tp of subharmonic modu¬ 
lation was investigated in two series of experiments 
on fiber lasers by Newell et al. [Newell et al., 1997], 
Celet et al. [Dangoisse et al, 1997] and one on a 
CO 2 laser by Chizhevsky et al. [1986]. The former 
concentrates on checking the predictions of the the¬ 
ory presented in the preceding paragraph, i.e. on the 
dynamics in the vicinity of the first period-doubling 
transition while Celet et al. aims at a global but 
less precise investigation of the effect of this phase 
p>, i.e. shift and/or change of nature of all bifur¬ 
cation points versus the phase of the subharmonic 
modulation. Chizhevsky et al. carried out several 
investigations on the role of phase in a driven sys¬ 
tem subjected to harmonic and subharmonic mod¬ 
ulations with special emphasis on the bistability in¬ 
duced by the second modulation in the presence of 
a swept parameter. In their study the nonadiabatic 
effects due to the sweeping of the control parameter 
strongly interefere with those due to subharmonic 
modulation. We concentrate here on experiments 
and theory in which these effects are made negli¬ 
gible as in [Newell et al., 1997; Dangoisse et al., 
1997], 

Both series of experiments have been carried 
out on Neodymium-doped fiber lasers pumped by 


diode lasers operating at 810 nm. Most experi¬ 
mental characteristics of the experimental set-ups 
are similar: fiber length 2.8 and 3.5 m, mirror re¬ 
flecting coefficient 0.95, pump power 1.93 and 2 
times above threshold respectively, making these 
two works quite complementary. 

By probing the 82 — <p control space, Newell 
et al. could show that the response of the laser is 
drastically changed if <p = 0 and minimized when 
ip = 7 t/2 . The bifurcation diagrams they obtained 
are displayed on Fig. 3. In both sets the bifurca¬ 
tion diagram in presence of modulation is compared 
with a reference one obtained in absence of addi¬ 
tional modulation (82 = 0). The existence of two 
branches when 82 7^ 0 clearly shows that no true PI 
orbit exists. The difference in the splitting of the 
branches of the bifurcation diagram [Figs. 3(a) and 
3(b)] in the vicinity of the period-doubling bifurca¬ 
tion is clear and strongly dependent on the phase <p 
of subharmonic modulation in accordance with the 
predictions. 

Examples of this effect together with the mod¬ 
ifications of the other bifurcations are shown in 
Figs. 4 and 5. They show bifurcation diagrams ob¬ 
tained by numerical simulations with the cw part of 
the pump Aq as the bifurcation parameter on the 
basis of a model of the fiber laser which proved very 
efficient in our previous studies of the dynamics of 
this kind of lasers [Derozier et al., 1992]. In the 
following, experimental and numerical results have 
been obtained with a modulation of the form 

A = A 0 + A\ cos u)\t + A 2 cos(w 2 t + V 3 ) 

for technical reasons. Extreme care has been taken 
to avoid any transient regime and slow passage ef¬ 
fects in all simulations and experiments. In typical 
simulations we compare the bifurcation diagrams 
obtained with increasing and decreasing values of 
the control parameter to reveal the possible gener¬ 
alized bistability effects due to the coexistence of 
attractors, especially in the vicinity of boundary 
crises. These bistability effects are not studied in 
detail here and are partly discussed in a preceding 
paper [Dangoisse et al, 1997]. 

Figure 4 shows the influence of a change in the 
amplitude of the additional modulation in the case 
of the fourth subhamonic modulation {uj\ = 4 ^ 2 ). 
Figure 4(a) is a reference in presence of the main 
modulation only while Figs. 4(b) and 4(c) were ob¬ 
tained for relative modulation amplitudes of 5.10 * 
and 5.10 -3 respectively, the amplitude of the main 
modulation is 0.42 relative to the cw pump for all 
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Pump parameter 



Fig. 3. Experimental demonstration of the influence of the 
phase <p of second subharmonic modulation on the first 
period-doubling bifurcation for a fiber laser. Diagrams at 
the top are references in absence of modulation, and those at 
the bottom correspond to additional (a) P = 0, (b) ip = tt/2 
subharmonic modulation. 


Fig. 4. Bifurcation diagrams versus amplitude of the sub¬ 
harmonic modulation for a model of the fiber laser (a) ref¬ 
erence in absence of modulation r 2 = 0, (b) rn = 5.10 -4 , 
(c) V 2 — 5.10 -3 . 

these recordings. Note that in all these simula¬ 
tions subharmonic modulations are much weaker 
than the main driving field. Other parameters and 
details about the model are given in [Dangoisse 
et al, 1997]. The two branches of the original di¬ 
agram in the domain A = 1.80 — 1.855 split into 
four branches as expected in the case of n = 4 
subharmonic modulation and the 2 T — 4 T bifur¬ 
cation becomes an imperfect bifurcation. The split¬ 
ting of the 4T branches obviously increases with r 2 , 
it is barely visible for A = 1.80 and = 5.10 -4 
but becomes quite clear for = 5.10~ 3 . More 
interestingly, the whole bifurcation diagram shifts 
to lower pump values and chaos appears in a pa¬ 
rameter domain where the system displays periodic 
response in absence of the second modulation and 
vice-versa. For instance, the one-tone modulated 
system is 4T periodic at Ao = 1.875 and become 
8 T periodic and chaotic for r 2 = 5.10 -4 and 5.10 -3 
respectively. 

The additional modulation also broadens the 
period-doubling cascade as demonstrated by the 
comparison of Figs. 4(a) and 4(c). More precisely, 
the cascade from the 4T-8T bifurcation to the 
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(a) 


(b) 


(c) 


Pump parameter 


Fig. 5. Influence of the phase of subharmonic modulation on 
the bifurcation diagram for a model of the fiber laser (a) ref¬ 
erence diagram in the absence of subharmonic modulation, 
(b) in the presence of in-phase (</j = 0), and (c) dephased 
(<p = 27r/9) fourth subharmonic modulation. Diagrams cor¬ 
respond to increasing pump parameters. 


transition point to chaos extends on a range of A 
values which doubles in the presence of the second 
modulation, from A A ~ 0.06 in absence of modu¬ 
lation [Fig. 4(a)] to AA ~ 0.12 in the conditions of 
Fig. 4(c). As can be seen on Figs. 4(b) and 4(c), the 
second modulation also induces periodic windows in 
the parameter domain where the laser behavior was 
chaotic when r-i = 0. For instance this is the case 
for the simulations reported on Fig. 4(c) in the do¬ 
main A = 1.891 — 1.894 where the original system 
is chaotic but the perturbed one is periodic. 

However the main information which can be 
extracted from these simulations is that both the 
shift of bifurcation points and the position of 
periodic windows strongly depend on the phase of 
the subharmonic with respect to the main modula¬ 
tion as shown on Fig. 5. This figure compares bifur¬ 
cation diagrams obtained in absence of subharmonic 
modulation and with two different phases {<p = 0 
and 27 t/ 9) for this modulation with respect to the 
main one. The reference diagram [Fig. 5(a)] dis¬ 
plays a period-doubling cascade transition to chaos 


between A = 1.76 and 1.90 in absence of subhar¬ 
monic modulation. Addition of a weak fourth sub¬ 
harmonic modulation leads to smaller A values for 
the bifurcation points but the magnitude of the shift 
towards lower A values depends on <p. Here the 
width of the chaotic region is reduced. Periodic 
windows and associated period-doubling cascades 
appear in the domain A ~ 1.88 -1.90 where the un¬ 
perturbed laser was chaotic. This definitely shows 
that the phase (f may be used to control the dynam¬ 
ics. This control is very efficient in two respects: 
(i) the subharmonic modulation may be much 
weaker than the main one — typically 10 2 to 10 3 
times — and (ii) it allows control of almost all kinds 
of regimes, from periodic to chaotic and vice versa. 

The diagrams presented in Figs. 5(b) and 5(c) 
also show that subharmonic modulation induces or 
shifts crises. The most obvious effect is the emer¬ 
gence of the periodic windows discussed in the pre¬ 
ceding paragraph since these windows follow the 
disappearance of chaos in the boundary crises. In 
addition to that, more subtle effects show off at 
the end of the chaotic domain. In the absence of 
perturbation, the corresponding chaotic regime ex¬ 
tends up to A = 1.910 where a period 2 attractor 
appears. The position of this crisis is perturbed by 
the subharmonic pertubation and here the chaotic 
domain is extended to larger A values. The period 4 
regime which emerges at the crisis is simply due to 
the fact that in presence of n = 4 subharmonic mod¬ 
ulation, a 4T cycle replaces the 2 T attractor which 
was observed in this domain, as expected for low 
modulation amplitudes. 

All the above bifurcation diagrams have been 
obtained with increasing values of the control pa¬ 
rameter. Slightly different ones are sometimes ob¬ 
tained with decreasing A values, even when the slow 
passage effects are avoided because of the coexis¬ 
tence of attractors for the same values of the pa¬ 
rameters. The corresponding bistability effects are 
particularly visible near the boundary crises as usu¬ 
ally observed, e.g. in modulated loss CO 2 lasers. 
As expected the limits of these domains are re¬ 
gions of high sensitivity to additional modulation 
because the stability of the attractor which will dis¬ 
appear in the crisis is very low near these points. 
Note that crises also appear for lower A values and 
the scattered points obtained near A = 1.76 corre¬ 
spond to other periodic attractors which coexist in 
this parameter domain with the main one studied 
in this work. There is a small range of bistability 
between these two attractors as demonstrated by 
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bifurcation diagrams obtained with reverse (de¬ 
creasing A) sweep [Dangoisse et al, 1997]. 

Experiments carried out on a modulated fiber 
laser confirm the predictions of numerical simula¬ 
tions. Bifurcation diagrams have been recorded 
with the pump power as the bifurcation parame¬ 
ter and subharmonic ranks 2 and 4. In accordance 
with numerical simulations, they show that there is 
no strong dependence of the results on the subhar¬ 
monic rank u> 2 /uq. Experiments also confirm the 
predictions of numerical simulations in the follow¬ 
ing points (i) the additional modulation shifts and 
magnifies the whole bifurcation diagram, (ii) the 
phase of subharmonic modulation allows the control 
of dynamical regimes (iii) crises may be induced or 
suppressed by acting on the phase of subharmonic 
modulation. These are very common features which 
have been observed for a wide variety of parame¬ 
ters, some examples of which are given in [Dangoisse 
et al ., 1997] and will not be recalled here. 

Control of chaos has also been observed in 
presence of transition to chaos via quasiperiodic¬ 
ity [Berge et al, 1984], Period-doubling cascades 
are most easily observed for modulation frequen¬ 
cies equal to about half the relaxation frequency of 
the. free running laser. Quasiperiodicity is obtained 
in the fiber laser when the external modulation fre¬ 
quency is detuned from this parametric resonance 
condition. We have considered the effect of subhar¬ 
monic perturbation in these conditions. Contrary 
to the previous paragraphs, the effects are demon¬ 
strated here for fixed values of the laser parameters 
in the vicinity of quasiperiodic regimes. Dynamic 
regimes are identified by observation of the Poincare 
section of the experimental attractor as indicated 
in [Derozier et al, 1992], Quasiperiodic regimes are 
identified by a closed curve in the Poincare section, 
periodic ones by a discrete set of points and chaos 
by a strange attractor. As for the period-doubling 
cascade, the dynamics of the laser is altered by sub¬ 
harmonic modulation and it is possible to drive the 
laser in different regimes depending on the phase of 
this modulation. For instance in the conditions of 
Fig. 6 where the unperturbed (ri = 0.438, 7-2 = 0) 
modulated laser is quasiperiodic [Fig. 6(a)], it be¬ 
comes 8 T periodic with in phase (p = 0) modula¬ 
tion [Fig. 6(b)] while dephased (</? = 7 t/ 18) second 
subharmonic modulation with V 2 = 0.02 induces a 
chaotic regime [Fig. 6(c)]. We have also observed 
situations in which an initially chaotic laser is con¬ 
trolled to periodic regimes by acting on the phase of 
a weak subharmonic modulation [Dangoisse et al, 
1997], 



4 ( a -u-) 

Fig. 6. Control of chaos in the case of quasiperiodicity route 
to chaos in a fiber laser (a) reference quasiperiodic state, 
(b) cp = 0 second subharmonic modulation induces period 8 
regime, (c) <p = tt/ 18 second subharmonic modulation in¬ 
duces a chaotic regime. 


Experiments carried out on both period¬ 
doubling cascades and quasiperiodicity transitions 
to chaos have shown that these scenarios are in gen¬ 
eral globally preserved by the addition of a subhar¬ 
monic modulation but that the dynamical regimes 
observed for given parameters may be qualitatively 
different. Moreover there may be some dramatic 
changes in the chaotic regimes such as the creation 
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of new attractors. All these effects are highly de¬ 
pendent on the phase of subharmonic perturbation 
and allow the control of the dynamics of nonlinear 
systems with a very small action like a simple phase 
shift of a modulation acting on this system. This 
control is very versatile in the sense that in well- 
selected parameter regions, it is possible to choose 
from a wide variety of regimes by just acting on the 
phase of the subharmonic modulation. The main 
drawback of this approach is that it applies only to 
externally nonautonomous systems. 

6. Conclusion 

The introduction of the concept of control of chaos 
was made on a very general model based on the lo¬ 
cal dynamics, i.e. on a linear approximation of the 
dynamics in the vicinity of the unstable period or¬ 
bits. Here we have explored a different path which 
relies on a fully analytic description of the dynamics 
and showed that it is possible to obtain quantitative 
information on the control range, i.e. the parameter 
domain in which this control is efficient. These re¬ 
sults have been supported by both numerical simu¬ 
lations and experiments on different lasers. Numer¬ 
ical simulations have allowed for the consideration 
of the bifurcations beyond the first period-doubling 
transition. They have demonstrated that control of 
chaos through the phase of a subharmonic modula¬ 
tion is extremely efficient and applies to various bi¬ 
furcations including those of the quasiperiodic route 
to chaos. Experiments carried out on several lasers 
(YAG, fiber, CO 2 ) indicate that these methods are 
very general and consequently that this approach 
may be extended to other dynamical systems. In 
addition, the present work also provides an alter¬ 
native explanation to the results obtained for the 
control of chaos by subharmonic modulation with 
extremely small detunings (near resonant perturba¬ 
tion) [Vohra et al, 1995] since such detunings may 
be considered as a slowly drifting phase and extreme 
care should be taken to separate e.g. detunings as 
small as 0.1 Hz with a basic resonance frequency 
of 10 kHz from the effects discussed in the present 
paper. 
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This paper presents two control schemes for the chaotic dynamics of a CO 2 laser with feedback 
which can be applied after the recognition of a leading frequency of the motion in the power 
spectrum. The first one is realized by means of a selective feedback loop which rejects all the 
frequency components except that of the leading cycle to be stabilized. The second one consists 
in a resonant sinusoidal modulation of the control parameter. 


1. Introduction 

Lasers represent reliable systems to investigate non¬ 
linear phenomena such as chaotic dynamics, both 
from an experimental and a theoretical point of 
view. In fact, it was demonstrated by Haken [1975] 
that the Maxwell-Bloch equations for the dynam¬ 
ics of a two-level laser are equivalent to the Lorenz 
model describing convective turbulence in fluid 
dynamics [Lorenz, 1963]. A recent aspect of non¬ 
linear dynamics studies concerns the possibility of 
directing chaotic dynamics to periodic orbits or 
steady states by applying small external perturba¬ 
tions. Ott, Grebogi and Yorke [1990] proposed a 
general method to stabilize a given periodic orbit 
embedded in a chaotic attractor by means of small 
perturbations to a system parameter; such pertur¬ 
bations are proportional to the deviation of the sys¬ 
tem from the unstable fixed point, and this implies 
the presence of a suitable feedback loop (for re¬ 
views of the method and applications see [Shimbrot 
et al., 1993; Ditto et al., 1995; Petrov et al., 1993]). 
In the field of laser physics, the “occasional pro¬ 
portional feedback” method (OPF) [Hunt, 1991], 
derived from the OGY scheme, has been applied 
to stabilize periodic orbits and steady states of a 
chaotic multimode Nd:YAG laser with a nonlin¬ 
ear intracavity KTP (potassium titanyl phospate) 
crystal [Roy et al., 1992; Gills et al, 1992; Colet 


et al., 1994; Carr & Schwartz, 1995]. In such a 
case the stability regime has been extended over one 
order of magnitude with respect to the unperturbed 
system. The OPF has been also applied to con¬ 
trol chaotic frequency emission in lead-salt diode 
lasers [Chin et al, 1996]. Another variation of the 
OGY method, known as “minimal expected devia¬ 
tion” (MED), has been applied to stabilize periodic 
orbits of a NMR (nuclear magnetic resonance) laser 
[Reyl et al, 1993]. Bielawsky et al. [1993] used a 
feedback scheme proportional to the derivative of 
the laser intensity to stabilize the chaotic regime of 
a Nd 3+ doped fiber laser pumped by a single mode 
laser diode. 

At the same time, other experimental [Meucci 
et al, 1994; Liu et al, 1994, 1995; Chizhevsky & 
Glorieux, 1995; Chizhevsky et al, 1997] and theo¬ 
retical [Liu & Rios Leite, 1994; Colet & Braiman, 
1996; Vilaseca et al, 1996] works on lasers deal 
with different control schemes not based on feed¬ 
back loops, but still consisting of suitable small 
modulations of control parameters. 

An alternative strategy to achieve control of 
chaos is based on modifications of state variables 
instead of control parameters. This technique has 
been first introduced by Pyragas [1992], which pro¬ 
posed a correction signal proportional to the dif¬ 
ference between the values of a given variable at 
different times. The delay time is selected equal to 
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Fig. 1. Experimental setup. G: diffraction grating. LT: laser tube. EOM: electro-optic modulator. M: outcoupling mirror. 
D: HgCdTe detector. P: preamplifier. A: differential high-voltage amplifier amplifer. B: bias input. LOG: logarithmic 
converter. F: washout filter. The dotted lines represents the two control schemes: C \: feedback loop control (subharmonic 
filtering). C 2 : open loop control (resonant modulation). 


the period of the unstable orbit to be stabilized. 
Bielawski et al. [1994] have successfully applied this 
method to a modulated CO 2 laser. Simmendinger 
and Hess [1996] used the Pyragas scheme in a semi¬ 
conductor laser with optical feedback. In this work, 
we describe experimental testing of two methods for 
controlling chaotic dynamics in a CO 2 laser with 
electro-optic feedback, which is an autonomous sys¬ 
tem providing low-dimensional chaos. The first 
method involves a frequency domain approach by 
means of a filtering feedback loop where the only 
admitted frequency is that of the orbit to be sta¬ 
bilized [Genesio et al., 1993; Meucci et al., 1996; 
Ciofini et al., 1997]. The second method is based 
on the introduction of a small modulation of the 
control parameter at a frequency suggested by the 
power spectrum. 

The paper is organized as follows: In the sec¬ 
ond section, after a description of the experimental 
set up, we present the results obtained with the 
two methods. The third section is devoted to a 
presentation of the model and numerical simula¬ 
tions. Finally, we summarize the main results in the 
conclusions. 

2. Experimental Results 

The experimental setup is shown in Fig. 1. It con¬ 
sists of a single mode CO 2 laser with a feedback on 
the cavity losses, realized by means of an intracavity 
electro-optic modulator, driven by a signal propor¬ 
tional to the output intensity. The bias voltage B 
provided by the high-voltage amplifier is the control 
parameter of the system. The feedback voltage V 
obeys the following equation: 

V = A v - B + TGj) (1) 


where (5 = 300 kHz is the damping rate of the 
feedback loop, I is the adimensional laser inten¬ 
sity and R = 6.6 * 10~ 10 is the total gain of 
the feedback loop. The term al (a = 1.2 * 
10 -13 ) accounts for the nonlinearity of the detection 
apparatus, which consists of a HgCdTe detector. 
In this configuration, the laser undergoes a direct 
Hopf transition from a stationary to an oscillating 
regime beyond a critical bias value. For appropriate 
values of pump and feedback loop gain, the oscil¬ 
lating regime becomes chaotic through a sequence 
of subharmonic bifurcations. A three-dimensional 
reconstruction of the chaotic attractor, obtained for 
B = 360 V, is shown in Fig. 2. This reconstruc¬ 
tion has been performed by an embedding tech¬ 
nique of the laser intensity signal, using a delay time 
r = 3 p s. The delay time value has been chosen to 



Fig. 2. Three-dimensional reconstruction of the chaotic 
attractor for B = 360 V, r = 3 /is. 
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closely reproduce the projection of the attractor in 
the two-dimensional phase space I-V [Meucci et al ., 
1997]. We adopted a 3-D representation in Fig. 2 
instead of a 2-D representation because it provides 
better insight into the stretching and folding mech¬ 
anism which is necessary to maintain the chaotic or¬ 
bit within a finite volume of phase space. Another 
useful method to represent the dynamical proper¬ 
ties of a system is provided by Poincare sections. In 
Fig. 3 we show the Poincare map obtained by plot¬ 
ting an intensity maximum versus the previous one. 
The chaotic nature, namely the loss of information 
as time or iteration number increase, is related to 
the noninvertibility of the map 7(n + 1) = f(I(n)). 
In fact, due to this feature, it is always possi¬ 
ble to estimate J(n + 1) from I(n) (e.g. fitting 
the map of Fig. 3 with an appropriate polynomial 
expression), but there is ambiguity to retrieve I(n ) 
from I{n + 1). The fractal character of the attrac¬ 
tor is confirmed by an evaluation of the correlation 
dimension D 2 = 2.10 ± 0.04 [Grassberger & 

Procaccia, 1983]. Finally, an important feature re¬ 
lated to our control schemes can be extracted from 
the chaotic power spectrum of the laser intensity 
(Fig. 4), that is, the peak at around /o = 22 KHz, 
which is the remnant of the Hopf bifurcation fre¬ 
quency. 

The first control method has been realized by 
means of (Fig. 1). The transfer function of the 



Fig. 3. Poincare map, corresponding to the chaotic attrac¬ 
tor of Fig. 2, obtained by plotting an intensity maximum 
versus the previous one. 


filter, shown in Fig. 5, presents two zeroes: one cor¬ 
responding to the peak of the unperturbed chaotic 
spectrum /o, and the other at zero frequency. The 
maximum of the transfer function occurs at /o/2, 



Fig. 4. Power spectrum corresponding to the chaotic attrac¬ 
tor of Fig. 2. 



Fig. 5. Transfer function of the filter (a) Amplitude 
respones, (b) Phase response. 











Fig. 8. Modulation 
attractor (yellow). 



Fig. 9. Poincare maps corresponding 
period-4). 
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where the phase response crosses the zero. The ef¬ 
fect of these characteristics, when the filter is in¬ 
serted in the negative feedback loop of the laser 
intensity, is to stabilize the orbit corresponding to 
the only frequency component which is not fed back 
as a correction signal, that is the period-1 orbit 
with frequency /o embedded in the chaotic attrac¬ 
tor (Fig. 6). In this case the relative perturbation 
amplitude is 7%. If the controller is modified by 
inserting a logarithmic amplifier to drive the filter, 
its performance increases, providing stabilization of 
the period-1 orbit with smaller values of the relative 
perturbation. 

The Poincare maps of the unperturbed chaotic 
attractor and of the period-1 stabilized orbit are 
shown superimposed in Fig. 7. This feedback 
method also allows stabilization of the period-2 
orbit (Fig. 7). In this case however, the pertur¬ 
bation applied to the system is higher (12%) and, 
as a consequence, one of the points of the map does 
not lie exactly on the unperturbed map. Indeed, 
to optimize the control for the period-2 orbit (or 
for any other higher order subharmonics), one has 
to prepare a new filter with a zero at /o/2 and a 
maximum at /o/4 [Ciofini et al., 1997]. 

The second control method we have imple¬ 
mented (Fig. 1) consists of introducing a small sinu¬ 
soidal modulation of the control parameter B(t) = 
P * (1 + m * sin (u>t)) at a frequency close to that of 
the peak in the chaotic spectrum. Figure 8 
shows, superimposed on the chaotic trajectory, the 
stabilized period-4, period-2 and period-1 orbits 
obtained by increasing the relative perturbation 
amplitude m (5%, 8% and 22%, respectively) at 
a frequency of 21.6 kHz where the system presents 
the maximum sensitivity. Note that stabilization 
of the period-1 orbit requires a large amount of 
modulation. As a consequence, the orbit is af¬ 
fected by large harmonic distortions, and it is not 
embedded in the chaotic attractor. It is impor¬ 
tant to observe that the fundamental frequencies 
of the stabilized attractors are locked to the per¬ 
turbation frequency in agreement with the theory 
of periodic perturbations [Khalil, 1992]. The re¬ 
turn maps of the stabilized cycles are shown in 
Fig. 9 superimposed on the unperturbed chaotic 
map. The differences between the chaotic maps 
of Figs. 7 and 9, which however maintain the 
same general features, are due to the slightly 
different values of B for which they have been 
recorded. The differences in the performances of 
the two control methods are summarized in Table 1. 


Table 1. Relative perturbation amplitudes (exper¬ 
imental values) for the stabilization of period-1, 
period-2 and period-4 unstable orbits with the two dif¬ 
ferent control schemes. The selective filter, optimized 
for period-1 stabilization, does not allow stabilization 
of the period-4 cycle. 



Period-1 

Period-2 

Period-4 

Selective filter 

7% 

12% 

— 

Parametric mode. 

22% 

8% 

5% 


3. The Model 

The behavior of a single mode CO 2 laser is quantita¬ 
tively described by a set of five differential equations 
which, besides the radiative coupling between the 
resonant molecular transition (population inversion 
N ‘2 — ) and the field intensity, account also for the 

collisional transfer from the manifolds of the other 
rotational levels (population inversion M 2 — Mi). 
Considering the presence of the feedback [Eq. (1)], 
and a suitable rescaling of the variables, the model, 
also described in [Meucci et al., 1997], is: 

w = ko(x 2 — 1 — ki sin 2 (a^)) 
x 2 = -Tix 2 - 2 k 0 x 2 e w + 7x3 + x 4 + P 0 
x 3 = -rix 3 - x 5 + 7x2 + P 0 

( 2 ) 

X4 = —F2X4 - 7X5 + zx 2 + zP 0 
X5 = -r 2 x 5 - zx 3 + 7X4 + zPo 
x 6 = -/3x 6 + j3B 0 - f3f(e w ) 

where w = log(xi) (xi is the rescaled intensity), x 2 
is proportional to the population difference N 2 — Ni, 
X 3 to N 2 + Ni, X 4 to M 2 — Mi, X 5 to M 2 + Mi 
and X 6 to the feedback voltage V. The nonlin¬ 
earity of the detector is contained in the func¬ 
tion /(x), while B 0 represents the rescaled con¬ 
trol parameter. The other parameter values are: 
k 0 - 28.57, h = 4.56, r x = 10.06, T 2 = 1.06, 7 = 
0.05, P 0 = 1.6 * 10 ~ 2 and p = 0.43. 

In the frequency domain, Eqs. (2) can be rep¬ 
resented in the Lur’e form, shown in Fig. 10. L(s), 
where s — iuj, represents the transfer function of 
a linear dynamical block, corresponding to the sec¬ 
ond, third, fourth and fifth equations. The first 
and the last equations of ( 2 ) represent the feed¬ 
back to the linear block. Note that the variable 
X 2 is not an accessible quantity in the experiment, 
while the variable w, neglecting the nonlinearity of 
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the detection process, can be obtained after a loga¬ 
rithmic amplification of the laser intensity. In this 
schematization, the total gain of the feedback loop 
R = 133.9 has been split in two parts; r = 0.1339 is 
the gain associated to the optical detector (in series 
with the preamplifier) and R/r = 1000 is the gain 
of the high-voltage differential amplifier. 

The two control schemes are shown in Fig. 10. 
The first one ( C\ ) consists in the cascade of a log¬ 
arithmic amplifier and a linear selective filter with 
transfer function C(s). C(s ) fulfils the requirements 
given in the previous section, i.e. presents two ze¬ 
roes at / = 0 and f = fo, and a maximum at /o/2. 
The analytical expression of C(s ) is {loq = 2nfo): 

r <~\- ks(s 2 + oi 0 2 ) 

“ 7-2\- 

f s 2 -Ka>os + j (s + n) 

The parameter values are: cjo = 0.2016, k = 3.5, 
C = 0.7 and /x = 0.8. 

The second controller ( C 2 ) is simply a modu¬ 
lation signal applied to the summing point, which 
results in a sinusoidal modulation of the bias 
term Bo- 

Figure 11 shows a three-dimensional recon¬ 
struction of the chaotic attractor of the unper¬ 
turbed system obtained from numerical integration 
of Eqs. (2) (j 3 = 223 V), with the superposition of 
the period-1 orbit stabilized with the selective filter¬ 
ing feedback C\. In Fig. 12 we show the Poincare 
maps corresponding to the unperturbed attractor 
and the period-1 (k = 3.5) and period-2 (k = 2.9) 
orbits. Finally, Fig. 13 shows the results obtained 
by applying the resonant modulation C 2 to the bias 
voltage. In this case, we observe that the stabi¬ 
lized period-2 and period-4 orbits also visit regions 
external to that of the chaotic attractor. This fact 
indicates some limitations of our model, which is 
able to satisfactorily reproduce the experimental re¬ 
sults in the case of the filtering feedback control, but 
provides only a qualitative representation of the sta¬ 
bilized orbit in the case of the resonant modulation 
control. 

4. Conclusions 

In this paper we have shown that the chaotic 
dynamics of a CO 2 laser with feedback can be 
stabilized by using two different strategies, based 
on selective filtering and parametric modulation, 
respectively. From the results presented we summa¬ 


rize the main differences between the two proposed 
control schemes. On the one hand, the first method 
works correctly only for stabilization of period-1 cy¬ 
cle, and thus it seems less flexible than the second 
one, which provides stabilization of different peri¬ 
odic orbits only by changing the perturbation am¬ 
plitude. On the other hand, the relative perturba¬ 
tion amplitude introduced by the second method 
to stabilize the period-1 orbit is very large, so that 
the controlled trajectory is quite different from that 
embedded in the chaotic dynamics. 

As a final remark, it is important to com¬ 
pare the filtering feedback method with the time- 
delayed autosynchronization method proposed by 
Pyragas [1992]. In the frequency domain, the Pyra- 
gas scheme corresponds to performing a high-pass 
filtering without affecting the frequency component 
fo and its harmonics (all zeroes in the transfer 
function), thus ensuring that the stabilized orbit is 
exactly that embedded in the chaotic attractor. 
Anyway, at variance with our method, the presence 
of high frequency feedback makes the system sen¬ 
sitive to noise and reduces its robustness [Meucci 
et al., 1997]. 
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Based on our earlier theoretical work on tracking periodic patterns into spatiotemporal chaotic 
regimes in spatially extended systems, we present in this article a case study of a high-aspect- 
ratio three-level laser. Following a detailed investigation into travelling wave solutions of the 
laser system and their stability conditions, we discuss ways of stabilising these solutions by local 
feedback algorithms. In numerical simulations, we choose Pyragus continuous delayed feedback 
algorithm to be the local feedback form and demonstrate that stable travelling wave solutions 
of the spatially extended three-level laser can be greatly extended to unstable regions in the 
presence of this feedback. 


1. Introduction 

Spatiotemporal chaos occurs when different types 
of motion, excited in local regions in an extended 
system, interact to destroy the spatial coherence 
of the system concurrent with the onset of tempo¬ 
ral chaos. This phenomenon in continuous phys¬ 
ical systems is described by partial differential 
equations. While the transition from coherence to 
spatiotemporal chaos has yet to be characterised 
by global quantitative laws, certain normal mode 
equations have shown that such a chaotic state in 
a spatiotemporal context underlies different unsta¬ 
ble coherent structures, which are sensitive to small 
perturbations. The possibility of controlling spa¬ 
tiotemporal dynamics by these perturbations has 
recently inspired considerable theoretical and ex¬ 
perimental effort in many branches of nonlinear 
science [Sepulchre & Babloyantz, 1993; Hu & Qu, 
1994; Auerbach, 1994; Qin et al., 1994, Aranson 
et al, 1994; Johnson et al, 1995; Lourenco et al., 
1995; Petrov et al, 1995; Poon & Grebogi, 1995; 


Hagberg et al, 1996; Lu et al, 1996; Bleich & 
Socolar, 1996; Martin et al, 1996; Battogtokh & 
Mikhailov, 1996]. This offers an opportunity to 
stabilise, select and manipulate these systems for 
applications. 

To control a spatially extended dynamical sys¬ 
tem, a feedback with spatial coupling is often ap¬ 
plied to the system so that a desired state can be 
stabilised. Systems investigated are both one- and 
two-dimensional, continuous or discrete, in various 
disciplines. Based on these developments, control 
methodology has been further explored where em¬ 
phasis has been given to algorithms which are sim¬ 
ple and more readily implemented in experiments. 
In our recent work, it has been established that un¬ 
stable periodic solutions of spatially extended two- 
dimensional systems, when weakly perturbed, can 
be stabilised by local feedback without spatial cou¬ 
pling [Lu et al, 1997]. This is due to the fact 
that the evolution equations for weak perturba¬ 
tions to these periodic solutions can be reduced 
to ordinary differential equations, because of the 
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decoupling of the wave numbers in the perturba¬ 
tions to a first order small signal approximation. 
As a result, stable periodic solutions can be tracked 
to unstable and spatiotemporal chaotic regimes by 
simply utilising feedback control algorithms estab¬ 
lished for controlling temporal chaos in spatially 
constrained systems. In this article, we provide a 
case study of a high-aspect-ratio three-level laser 
in which a detailed investigation into the stability 
of travelling wave solutions of the laser system and 
ways to stabilise these solutions by local feedback 
are discussed. In numerical simulations, we choose 
the well-known Pyragus continuous delayed feed¬ 
back algorithm to be the local feedback and demon¬ 
strate that stable travelling waves of the spatially 
extended laser system can be greatly extended to 
unstable regions in the presence of this feedback. 

2. Theory for the Tracking Procedure 

We first give a brief review of the theory for tracking 
procedure by using local feedback algorithms. We 
consider a general form of a distributed nonlinear 
optical system 

J t '=N(q,n) + iDVlq , (1) 

which admits a travelling wave solution of 

go = C' e < < k - r - wt >, (2) 

where g is the complex amplitude of the electromag¬ 
netic field, N the nonlinear function describing the 
local field-material interaction, Vj_ the transverse 
Laplacian in two-dimensional space r(x, y ) and t 
the time, /i is the control parameter of the system 
and D a coefficient describing diffractive coupling 
in space, k and u> are the characteristic wave vec¬ 
tor and frequency of the travelling waves while C is 
the amplitude. The stability of the travelling solu¬ 
tion is determined by standard perturbation analy¬ 
sis which yields 

= iDV\5q + N'(q 0 , p)5q 

+ 2 N"(qo, fi)Sq 2 + 0(8q 3 ) , (3) 

where 5q is the perturbation strength and N' and 
N" are the first and second derivatives of N with 
respect to q at q = go- We note that terms of 5q* 
should appear in the above equation if N comprises 


explicitly the variable g*. In Fourier space, the per¬ 
turbation is of the form 

r+ oo 

Sq( r, t) = e *( k ' r ~ wt ) / 5q p e ip ' r dp , (4) 

J —OO 

where Sq p is the perturbation strength with wave 
vector p. By substituting Eq. (4) into Eq. (3), we 
derive the evolution equation of perturbation for the 
specific p 

= %u5q p -iD{k+p) 2 5q p +N' 8q p 

1 f + OO 

+ -N"e-^ / 6q p '5q p ^ p ,dp' + 0(5ql ). 

A J — OO 

( 5 ) 

We note that since Eq. (2) is a solution of Eq. (1), 
N'(qo, n), the Jacobian which along with k, p and 
u determines the linear stability of the solution, is 
independent of the transverse coordinate r. 

When the perturbation against the travelling 
wave solution is weak, the stability of this solution 
is determined by the linearisation of Eq. (5) which is 
a wave vector decoupled ordinary differential equa¬ 
tion. Control of such a solution, when unstable, can 
therefore be achieved by a feedback approach which 
is no more complicated than controlling an unstable 
periodic orbit in an ordinary differential equation; 
this is 

= iu8q p - iD (k + p ) 2 5q p + N'8q p + f p (t ), 

( 6 ) 

where f p is the feedback which depends only on the 
wave vector component p in the Fourier space. In 
r space the feedback is given as 

r+oo 

F(r, t) = e *( k ’ r-wt ) / f p (t)e p r dp (7) 

J — OO 

which does not involve feedback coupling in differ¬ 
ent spatial regions. Thus, the feedback required 
for this case is local in the transverse space, in the 
sense that the signal at a given point depends on 
only its behaviour at that particular point, and not 
on that in the neighbouring and distant regions. It 
follows that when a state q is close enough to a tar¬ 
geted travelling wave solution of a distributed sys¬ 
tem as given by Eq. (1), this state can be stabilised 
to the solution by simply adopting the control algo¬ 
rithms developed for controlling temporal chaos in 
ordinary differential equations. In order to utilise 
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the linearised approach to stabilise a traveling pat¬ 
tern, we may take advantage of the tracking proce¬ 
dure [Gills et al, 1992] which was first developed 
for extending a stable operating region of a low¬ 
dimensional laser system. In this procedure, the 
feedback to Eq. (1) is applied when the system is 
initially set in a region of a stable periodic pattern, 
and on varying a control parameter of the system, 
this stable solution can be extended to unstable and 
spatiotemporal chaotic regions in the presence of 
the feedback. Moreover, this tracking procedure 
can be adopted to manipulate a stabilised travelling 
wave pattern; the values of its frequency and wave 
vector can be altered, as can the orientation of the 
wave vector, by judiciously varying these parame¬ 
ters providing that they constitute one solution of 
the travelling wave family. The above analysis can 
be readily extended to stabilising more complex pe¬ 
riodic patterns, such as square and hexagonal pat¬ 
terns, in a set of coupled equations. 

3. Three-level Laser Model 

We consider a coherently optically pumped three- 
level laser as a case for study. It has previously 
received considerable attention as a spatially con¬ 
strained dynamical system in the investigation of 
temporal chaos [Moloney et al, 1992; Forysiak 
et al, 1991]. When the transverse space is broad, 
self-diffraction of the inter-cavity field plays an 
essential role. With inclusion of this transverse 
coupling the Maxwell-Bloch equations describing 
resonant single mode lasing are generalised to 

dE 

— = <7E + gP 2 + iaV 2 ± E , 

dPi 

~ = ~Pi+ AN x - EP 3 , 

dP 2 

~^r = -P 2 + EN 2 -APZ, 

dP 3 ' ( § ) 

1 JL = -p 3 + e*p 1 + api, 

^ = -b(l + N 1 )-2A(P l +P 1 *) 

- (E*P 2 + PP 2 *), 

dNo 

—■ = -WV 2 - A(Pt + Pf) - 2{E*P 2 + PP 2 *), 

where E is the slowly varying electromagnetic field 
amplitude of the laser emission and Pi, P 2 and 
P 3 are the normalised off-diagonal density matrix 


elements of the polarisation. Due to the exis¬ 
tence of transverse diffraction, these variables are 
complex even for the case of resonant pumping 
and laser emission, showing the important role of 
phase in optical pattern-forming phenomena. The 
other variables are Ni and N 2 , the diagonal den¬ 
sity matrix elements describing the population dif¬ 
ferences between the common upper and lower 
levels of the pump and lasing transitions. Pa¬ 
rameters of the system are a, the cavity damp¬ 
ing constant, g, the unsaturated gain of the laser 
medium, and a, the diffractive coefficient which can 
be set to unity by rescaling the transverse coordi¬ 
nates (x and y). Control parameters are A, the 
external pump strength, and b, the ratio of en¬ 
ergy relaxation ( 7 ) to dipole dephasing (T) rates 
of the medium. The laser system possesses trans¬ 
lational symmetry and is invariant under the fol¬ 
lowing transformation (E, P\, P 2 , P 3 , N\, N 2 ) => 
(-E, Pi, — P 2 , — P 3 , N \, N 2 ), which give certain 
restrictions to the solutions of Eqs. ( 8 ). 

3.1. Travelling wave solutions 

Equations ( 8 ) admits a travelling wave solution 

E = E°e i( ^' r - Ut \ Pi = Pi°, 

P 2 = P 2 V( k ' r ~ wt ), P 3 = pO e -i(k-r-'M ), (g) 

Ni = Nl N 2 = N 2 , 

where k and u> are the wave vector and frequency of 
the travelling wave solution on the transverse plane 
(x, y ), and E°, Pf, P 2 °, P 3 , N° and N 2 are the 
amplitudes, of which E° can be assumed to be real. 
By using the transformation E° = x\, u) = x 2 , 
Pi = Xz + ix 4 , P 2 ° = £5 + ix 6 , P° = X7 + ix 8 , 
= x 9 , N 2 = xw and substituting Eqs. (9) into 
Eqs. ( 8 ), we derive the following coupled equations 
for x\, x 2 , Xg and xw for a nontrivial solution, 

(1 + x[)\gx 10 - a - x 2 {k 2 - x 2 )\ 

- x 2 [k 2 - (1 + a) x 2 ] - A 2 (cr + gx 9) = 0 , 

(1 + zf )\k 2 - (1 + a)x 2 ] + x 2 [gxio - a 

- x 2 (k 2 - x 2 )} + A 2 (k 2 - x 2 ) = 0 , 

2A 2 xg + (b- 2x\)x 10 (10) 

2 

+ -x\ [3cr + x 2 (k 2 - x 2 )\ — 0 , 
b+(b + 4 A 2 )x g - 4x 2 xio 

+ ^x 2 [3a -t- 2a - 2 (k 2 - x 2 )} = 0, 
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while other variables are given by 
£3 = Ax 9 — X\Xj , 


trial solution 


£4 = —0:1X8 ; 
a 

£5 = -xi, 

9 


x e = -(k 2 -x 2 )x 1 , v 11 ) 

9 

x 7 = -Jgl9 x 10 - <7 - x 2 (k 2 - x 2 )}, 

x 8 = ~[k 2 - {1 + a)x 2 \. 

Ag 

Figure 1 shows a set of travelling wave solutions in 
three-dimensional space; the laser amplitude is the 
function of the pump strength A and the wave num¬ 
ber k, with other parameters being held as g = 52, 
b = 0.4 and a = 1.3. The curve for k = 0 corre¬ 
sponds to homogeneous steady state solutions, the 
laser branch (E A 0) of which exists only in a lim¬ 
ited region of the pump strength due to the Rabi 
splitting effect of the laser gain profile. On increas¬ 
ing the value of k form a family of nontrival trav¬ 
elling wave solutions, which exist over an extended 
region of the pump strength but is limited only for 
small wave numbers. 


3.2. Stability of travelling 
wave solutions 

In general, the stability of the traveling solutions 
given by Eqs. (9) through Eqs. (10) and (11) de¬ 
pends on the control parameters of the laser system. 
To determine the stability conditions of these trav¬ 
elling waves, we perturb them with the following 



/ +00 1 

8E p e ip - r dp e i ( k ' r_a '*), 

-OO 

/ +OO 

8P\ p e ip ' T dp , 

-OO 

/*+oo "I 

P 2 = po+ / 8P 2p e ipr dp e i ( k ' r ~ tJi ), 

Z • ( 12 ) 

P 3 = P 3 0 + / 8P 3p e ipT dp e-’( k ' r - ut ), 

J —OO 

/ +OO 

8N\ p e ip ' T dp , 

-OO 

/ +OO 

8N 2p e ipr dp , 

-OO 

where the second terms on the right-hand side are 
perturbations integrated over all wave numbers. By 
substituting Eqs. (12) into Eqs. ( 8 ) and considering 
only to the first order small terms, we derive the 
following linearised evolution equation for the per¬ 
turbation strength 

= iu5E p - i (k + p ) 2 5E p - cr8E p + g8P 2p , 


l -E = —5P\p + A5Ni p - E°8P Zp - P°5E p , 


d5P 2p 


ddP ?J p 


= -b8Ni p - 2 A(SPi p + 5P{ V ) - ( E°8P 2p 
+ P 2 8E P + E°5P* p + P°*SE p ), 

= -b8N 2p - A(8P lp + 8P { p ) - 2 (E°5P 2p 
+ P°8E p + E°5P; p + P°*6E p ), 

where (k + p) 2 is equal to (k+p) 2 for the two vec¬ 
tors being parallel, corresponding to the Eckhaus 
instability, and to k 2 +p 2 for the vectors being per¬ 
pendicular, a case of Zig-zag instability. The eigen¬ 
value equations for X p are obtained by substituting 
the relation d8( )/dt = A p 8( ) into Eqs. (13), where 


= iu8P 2p — 8P 2p + E° 8 N 2p 


+ N°8E P - A8P 3 * p , 


= —iu8P 3p — 8P;i P + E° 8 P\ 
+ p?se; + A8P 2p , 


Fig. 1. Travelling wave solutions, B as a function of the 
pump strength A and wave number k for parameters g = 
52.0, b = 0.4 and a = 1.3. 
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Fig. 2. (A, k) two-dimensional stability diagram for trav¬ 

elling wave solutions for parameters g = 52.0, b — 0.4 and 
(T = 1.3. The region surrounded by the solid curves comprises 
nontrival solutions in which domains marked by El and ZI 
correspond to Eckhaus and Zig-zag instabilites. N and S 
stands for regions of trival solutions and stable solution re¬ 
spectively. 

8( ) stands for the above perturbations. The travel¬ 
ling wave solutions are linearly stable if Re(A) < 0 
and unstable if Re(A) > 0 . Figure 2 shows a bi¬ 
furcation diagram of the steady state and travelling 
wave solutions in (A k ) two-dimensional space, for 
parameters g = 52.0, b = 0.4 and a = 1.3. The re¬ 
gion surrounded by the solid curves is that of non¬ 
trival travelling solutions in which domains marked 
by El correspond to the Eckhaus instabilites and 
the area under ZI is Zig-zag unstable. N stands for 
regions of trival solutions. 

4. Tracking of Travelling 
Wave Solutions 

In this section, we apply the general theory in Sec. 2 
to the laser equations ( 8 ) to track the travelling 
wave solutions from a stable region into an unstable 
domain. As we discussed earlier, tracking of such 
unstable solutions can be realised by simply adopt¬ 
ing the control approaches developed for controlling 
temporal chaos in ordinary differential equations. 
Well-known methods for control of temporal chaos, 
such as the Ott-Grebogi-Yorke algorithm (OGY) 
[Ott et al. , 1990], occasional proportional feedback 
(OPF) [Hunt, 1991] and continuous delayed feed¬ 
back (CDF) [Pyragas, 1992], can therefore be used. 
Here we take the continuous delayed feedback as an 



Pump Strength A 

Fig. 3. (a) Travelling wave (tw) and homogeneous steady 

state (hs) solutions as a function of the pump strength for 
g — 52.0, b = 0.4, o = 1.3 and k — 1.78. The solid and 
dash curves stand for stable and unstable solutions respec¬ 
tively, separated by the Hopf bifurcation (HB) points, which 
for the travelling wave is further referred to as the Eckhaus 
instability (ECK). (b) The corresponding wave numbers and 
frequencies of the travelling wave in trace(a), both as the 
function of the pump strength. 

example. For this case, the form of the feedback as 
discussed in Eq. ( 6 ) is f p ~ ( 8E p {t - to) - 8E p {t)). 
Using Eq. (7) the feedback in r space is given as 

F(r, t) = a{E( r, t-t 0 )~ E{ r, t)) , (14) 

which is now applied to the right-hand side of the 
first equation of Eqs. ( 8 ). In the above equation a is 
the proportionality constant defining the feedback 
strength and to the period in time of the travelling 
wave solution to be stabilised. 

We have simulated the tracking procedure by 
numerically integrating the three-level laser system 
in the presence of CDF control in the laser field 
equation. The numerical integration is based on 
the Fourier transformation method by using typi¬ 
cally 60 x 60 grid points in a square in the trans¬ 
verse plane and is checked by 120 x 120 grid point 
simulations. The width of the square is chosen to be 
32n/kA, where Fq is the spatial wave number cor¬ 
responding to the largest growth rate of the unsta¬ 
ble travelling wave solution. An example is shown 
in Fig. 3. The travelling wave solutions, the fre¬ 
quencies and wave numbers of which are given in 
trace (b), are stable only in a small pump region as 
depicted in trace (a). These solutions corresponds 
to the vertical line of k = 1.78 in Fig. 2 in which 
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the stable region is bounded by two points of the 
Eckhaus instability. To track this travelling wave 
branch, the feedback is first applied to each grid 
point, when the system is initially in the parame¬ 
ter region for stable travelling waves. The pump 
strength A is then increased from this region in 
small steps 8A — 0.1 to extend to the unstable re¬ 
gion. At each value of the pump strength, the initial 
conditions of the system are set to be the output of 
the system at "its previous value, whereas the delay 
time to in the feedback is chosen to be the period 
of the travelling wave at the current value of the 
pump strength. Numerical simulations have shown 
that the travelling wave can be stabilised along the 
dash curves of k and uj in Fig. 3(b), extending the 
periodic pattern into both the spiral and spatiotem- 
poral chaos regions. For each value of the pump 
strength, there is a threshold of the coefficient a for 
a successful tracking, the values of which are less 
than 1.0. In the process of tracking, the stabilised 
travelling waves are found to be robust against ex¬ 
ternal perturbations to the pump strength A. For 
a = 0.7 and different values of A, the system is 
shown to be stable when A is perturbed by differ¬ 
ent levels of white noise signal (peak value) up to 
20% of the signal. We further find that the feedback 
strength throughout the tracking process is small, 
demonstrating the noninvasive nature of the con¬ 
trol. In these simulations it is limited to less than 
2 % of the laser output (peak to peak) due to the 
limited temporal and spatial resolution in the sim¬ 
ulation scheme. We note that for stabilizing a state 
when it is far away from a targeted periodic pattern, 
the feedback by contrast requires both temporal and 
spatial coupling [Lu et al, 1996]. 

The tracking procedure has also been tested nu¬ 
merically on travelling wave solutions of different 
sets of wave numbers k and frequencies u. More¬ 
over, such a tracking procedure can be applied to 
track a travelling wave solution using other control 
parameters of the system, for instance, the wave 
number k, while the pump strength is held con¬ 
stant. For this case, the wave number of a stabilised 
travelling wave can vary while the corresponding 
frequency is judiciously chosen so that these two 
parameters satisfy solutions of the travelling wave 
family. In such a way, the spectrum of the peri¬ 
odicity of the stabilised pattern can be extended. 
In summary, we have shown that, as an exam¬ 
ple of a spatially extended dynamical system, the 
high-aspect-ratio three-level laser is well within our 
general theoretical description for the use of 


the tracking procedure to extend its coherent pat¬ 
tern formations. The advantage of this tracking 
procedure lies in the use of local feedback in a spa¬ 
tially extended system. Such a simple control ap¬ 
proach with no involvement of spatial coupling is 
more easily accessible to experiments. In optics, the 
CFD algorithm can be implemented using the inter¬ 
ference technogue of theoutput beam with delayed 
itself. 


Acknowledgments 

This work is supported by EPSRC (UK) Grants 
No. GR/J73285 and No. GR/K23768. 


References 

Aranson, I., Levine, H. k Tsimring, L. [1994] “Con¬ 
trolling spatiotemporal chaos,”' Phys. Rev. Lett. 72, 
2561-2564. 

Auerbach, D. [1994] “Controlling extended systems of 
chaotic elements,” Phys. Rev. Lett. 72, 1184-1187. 

Battogtokh, D. k Mikhailov, A. [1996] “Controlling tur¬ 
bulence in the complex Ginzburg-Landau equation,” 
Physica D90, 84-95. 

Bleich, M. B. k Socolar, J. E. S [1996] “Controlling 
spatiotemporal dynamics with time-delay feedback,” 
Phys. Rev. E54, R17-R20. 

Forysiak, W., Moloney, J. V. k Harrison, R. G. 
[1991] “Bifurcation of an optically pumped 3-level 
laser model,” Physica D53, 162-186, and references 
therein. 

Gills, Z., Iwata, C., Roy, R., Schwartz, I. B. k Triandaf, 
I. [1992] “Tracking unstable steady-state — Extend¬ 
ing the stability regime of a multimode laser system,” 
Phys. Rev. Lett. 69, 3169-3172 

Hagberg, A., Meron, E., Rubinstein, I. k Zaltzman, B. 
[1996] “Controlling domain patterns far from equilib¬ 
rium,” Phys. Rev. Lett. 76, 427-430. 

Hu, G. k Qu, K. [1994] “Controlling spatiotemporal 
chaos in map lattice systems,” Phys. Rev. Lett. 72, 
68-71. 

Hunt, E. R. [1991] “Stabilizing high-period orbits in a 
chaotic system — The diode resonator,” Phys. Rev. 
Lett. 67, 1953-1955. 

Johnson, G. A., Locher, M. k Hunt, E. R. [1995] “Sta¬ 
bilised spatiotemporal waves in a convectively unsta¬ 
ble open flow system — Coupled diode resonators,” 
Phys. Rev. E51, R1625-1628. 

Lourenco, C., Hougardy, M. k Babloyantz, A. [1995] 
“Control of low-dimensional spatiotemporal chaos in 
Fourier space,” Phys. Rev. E52, 1528-1532. 

Lu, W., Yu, D. k Harrison, R. G. [1996] “Control of pat¬ 
terns in spatiotemporal chaos in optics,” Phys. Rev. 
Lett. 76, 3316-3319. 








Instabilities and Tracking of Travelling Wave Patterns in a Three-Level Laser 1775 


Lu, W., Yu, D. k Harrison, R. G. [1997] “Tracking pe¬ 
riodic patterns into spatiotemporal chaotic regimes,” 
Phys. Rev. Lett. 78, 4375-4378. 

Martin, R., Scroggie, A. J., Oppo, G.-L. k Firth, W. J. 
[1996] “Stabilization, selection, and tracking of unsta¬ 
ble patterns by Fourier space techniques,” Phys. Rev. 
Lett. 77, 4007-4010. 

Moloney, J. V., Forysiak, W., Uppal, J. S. k Harrison, 
R. G. [1989] “Regular and chaotic dynamics of op¬ 
tically pumped molecular lasers,” Phys. Rev. A39, 
1277-1285. 

Ott, E., Grebogi, C. k Yorke, J. A. [1990] “Controlling 
chaos,” Phys. Rev. Lett. 64, 1196-1199. 

Petrov, V., Metens, S., Borckmans, P., Dewel, G. k 
Showalter, K. [1995] “Tracking unstable Turing pat¬ 


terns through mixed-mode spatiotemporal chaos,” 
Phys. Rev. Lett. 75, 2895-2898. 

Poon, L. k Grebogi, C. [1995] “Controlling complexity,” 
Phy. Rev. Lett. 75, 4023-4026. 

Pyragas, K. [1992] “Continuous control of chaos by self 
controlling feedback,” Phys. Lett. A170, 421-428. 

Qin, F., Wolf, E. E. k Chang, H. C. [1994] “Controlling 
spatiotemporal patterns on a catalytic wafer,” Phys. 
Rev. Lett. 72, 1459-1462. 

Sepulchre, J. A. k Babloyantz, A. [1993] “Controlling 
chaos in a network of oscillations,” Phys. Rev. E48, 
945-950. 







International Journal of Bifurcation and Chaos, Vol. 8, No. 9 (1998) 1777-1782 
© World Scientific Publishing Company 


EXPERIMENTAL SWITCHINGS IN BISTABILITY 
DOMAINS INDUCED BY RESONANT PERTURBATIONS 

V. N. CHIZHEVSKY*>t>*, R. VILASECA* and R. CORBALANt 
* Departament de Fisica i Enginyeria Nuclear, Universitat Politecnica de Catalunya, 

Colom 11, E-08222 Terrassa, Spain 
t Departament de Fisica, Universitat Autonoma de Barcelona, 

E-08193 Bellaterra, Barcelona, Spain 
*Stepanov Institute of Physics, Belarus Academy of Sciences, 

220072 Minsk, Belarus 

Received July 31, 1997; Revised October 28, 1997 


We show that a combination of a resonant perturbation which induces bistability in the system 
with the targeting technique based on the action of a short-lived pulse perturbation makes 
the nonfeedback control of nonlinear systems (not only in a chaotic state) more flexible and 
even competitive with OGY’s method in the sense of fast switching between controlled orbits 
belonging to coexisting attractors. For different initial states we present several experimental 
and numerical examples of such a type of control, which does not require feedback system. 


1. Introduction 

Controlling chaos and, more generally, controlling 
dynamics and complexity of nonlinear systems 
attracted much attention over the last years. 
For these purposes a diversity of approaches and 
methods has been proposed. Many of them have 
been successfully implemented experimentally in 
diverse systems ranging from physics, and optics to 
biology and chemistry and showed high efficiency 
and flexibility. The latter is especially true with 
regard to the method of chaos control proposed by 
Ott, Grebogi and Yorke [Ott et al., 1990; Ditto 
et al., 1995]. Among the different approaches there 
is the so-called nonfeedback method which is based 
on the adding of weak periodic resonant perturba¬ 
tions to the system under control [Lima k Petttini, 
1990; Braiman k Goldhirsch, 1991]. Typically, 
resonant perturbations at a subharmonic frequency, 
depending on their amplitude and phase and 
on the initial state, may produce stabilizing or 
destabilizing effects on the nonlinear system. As it 
has been experimentally shown, they may suppress 


chaos [Azevedo k Rezende, 1991; Fronzoni et al, 
1991; Ding et al., 1994; Meucci et al, 1994], 
suppress periodic orbits, induce a chaotic behavior 
if the system is initially in a periodic state, and, 
finally, may induce crisis of strange attractors 
[Chizhevsky k Corbalan, 1996]. Up till now, most 
attention was paid to the question of how resonant 
perturbations modify the dynamics of nonlinear 
systems. Recently, it has been experimentally and 
numerically shown with a loss-modulated CO 2 laser 
that resonant perturbations globally change the 
phase space of the system, splitting the primary 
attractor into two new ones [Chizhevsky et al, 
1997a]. It has been demonstrated that by changing 
the amplitude and/or the phase of resonant pertur¬ 
bations, one can control the overlapping of different 
dynamical regimes associated with these two new 
attractors. From the practical point of view, for 
example, in engineering, a new problem arises here: 
How to switch the motion in the system from 
one dynamical regime to another one belonging to 
these new attractors? Recently, the technique of 
targeting periodic orbits has been proposed and 
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developed experimentally [Chizhevsky &: Turovets, 
1994; Chizhevsky & Glorieux, 1995], which allows 
one (i) to switch the phase of a periodic motion, 

(ii) to reach unstable orbits created either at a 
period-doubling or saddle-node bifurcation, and 

(iii) to perform switching between stable periodic 
orbits belonging to coexisting attractors. It is based 
on the action of a large-amplitude pulse perturba¬ 
tion to the system that is equivalent to a change of 
initial conditions. There exists the optimal timing 
and the optimal amplitude of a pulse perturbation 
which allow one to perform switchings practically 
without transients [Chizhevsky et al., 1997b]. 

The main aim of this work is to show 
that a combination of the resonant perturbation, 
which induces bistability in the system, with the 
targeting technique makes the nonfeedback control 
of nonlinear systems (not only in a chaotic state) 
more flexible and even competitive with OGY’s 
method in the sense of fast switching between 
controlled orbits belonging to coexisting attractors. 
For different initial states we demonstrate several 
experimental and numerical examples of such a 
type of control. This technique does not require 
any feedback system and might be experimentally 
implemented in diverse nonlinear systems. 

2. Experimental Setup and Model 

The experiments were carried out on a CO 2 laser 
with two acousto-optic modulators inserted in the 
laser cavity as described earlier [Chizhevsky &; 
Corbalan, 1996]. Two electric signals were applied 
to the modulator providing the time-dependent 
cavity losses. The driving signal had the frequency 
/d = 100 kHz and the amplitude V&. The per¬ 
turbation signal had the frequency / p = /d /2 = 
50 kHz, the amplitude V], and the phase tp. Laser 
parameters such as relaxation-oscillation frequency 
/ro = 50 kHz and decay rate 7 ro = 200 kHz 
were estimated from the laser response due to 
short-lived loss perturbations (in the absence of 
the driving and perturbation signals) and were 
used then in numerical simulations. The laser 
responses were detected with a CdHgTe detec¬ 
tor and a digital oscilloscope coupled with a PC. 
The technique of stroboscopic data recording with 
a sampling period T (T = l//d) was used. 

Short-lived pulse loss perturbations were caused 
by absorption of the laser emission on nonequi¬ 
librium charge carriers which were excited by 
illuminating the semiconductor window of the laser 
tube by short pulses of a YAG-laser [Chizhevsky & 


Turovets, 1994; Chizhevsky Sz Glorieux, 1995]. 
Several clocks and devices were used in experiments 
to synchronize all control signals. 

Numerical simulations were performed using a 
simple two-level rate-equation laser model [Tredicce 
et al., 1986]: 


^ = t 1 {y-k)u, 

(1) 

^ = (yo - y)i ~ uy, 

( 2 ) 


where 

k = ko + fcj cos(27r ft) 

+ k p cos(TTft + ip) + ks(t - to) (3) 

Here u is proportional to the radiation density, y 
and j/o are the gain and the unsaturated gain in 
the active medium, respectively, r is half round- 
trip time of light in the cavity, 7 is the gain de¬ 
cay rate, k is the total cavity losses, ko is the con¬ 
stant part of the losses, k<± is the driving ampli¬ 
tude, k p is the perturbation amplitude, <p is the 
perturbation phase, k$(t — to) is the pulse loss per¬ 
turbation, to is the moment of time of the action 
of the pulse losses perturbation (ks(t — to) = 0 
for t < to and kg(t — to) = ks(exp(—a(t — to)) — 
exp {—(/3t — to))) for t > to ,where a = 3 x 10 7 
s -1 and (3 = 10 8 s -1 ). In order to compare 
the numerical results with the experimental ones 
the values /ro and 7 ro measured in experiments 
were used for finding some parameters appearing 
in Eqs. (1) and (2) using the following expressions: 
7 RO = 7 + uo (where u.q is a stationary value 
of u) and /ro = l/2n[(T~ 1 u 0 yo) ~ (7Ro/2) 2 ] 1/2 
which can be obtained from Eqs. (1) and (2) 
by a standard linearization procedure. Through¬ 
out our calculations the following fixed parameters 
were used: r = 3.5 x 10 -9 s, 7 = 1.978 x 10 5 , 
yo = 0.175, ko = 0.17303. The other parameters 
were varied in the simulations. Because in the 
numerical simulation and the experiments the con¬ 
trol parameter (k^ in the simulations and V ( \ in the 
experiments) is measured in different units, we put 
also to the figure captions the value of normalized 
control parameter p defined as p = k^/ki /2 — 
Vd/Vi/ 2 , where 2 and Vi / 2 are the values of the 
first period doubling threshold in the simulations 
and experiments, respectively. 

3. Results and Discussions 

The appearance of bistability induced by resonant 
perturbation is shown in Fig. 1 where numerical 
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Fig. 1. Numerical bifurcation diagrams of a CO 2 laser show¬ 
ing the appearance of the bistability, (a) The laser response 
without the resonant perturbation, (b) and (c) the laser re¬ 
sponses with the resonant perturbation. The parameters 
used in the computer simulation: /ro = 50 kHz, 7ro = 
200 kHz, /d = 80 kHz, fed varied between 1.24 x 10 -4 and 
4.4 x 10 -5 , (a) k p = 0, (b) and (c) k p = 9 x 10~ 6 , ip = 0. 


bifurcation diagrams of a CO 2 laser are presented. 
They were obtained with decreasing of the driving 
amplitude k<\ by a linear law. Shown here are 
maximal peaks of the laser intensity versus time. 
Figure 1(a) shows the laser response without the 
resonant perturbation. Figures 1(b) and 1(c) 
correspond to two laser responses which appeared 
instead of the primary solution in the presence of 
the periodic perturbation. It is seen [Fig. 1(b)] 
that one solution (the solution-1) is suppressed. As 
the driving amplitude decreases, the system reaches 
the first bifurcation point shifted by the resonant 
perturbation and the solution-1 becomes unstable. 
Then the system jumps to the solution-2. This 
solution is entirely presented in Fig. 1(c) (for more 
details see [Chizhevsky et al ., 1997a]). By com¬ 
paring Figs. 1(b) and 1(c) it is seen that for the 
given perturbation amplitude different dynamical 
regimes associated with solution-2 are overlapped 
with a 2 T stable regime belonging to the solution-1 


as the driving amplitude changes. The second ef¬ 
fect which one can see from Fig. 1 is that the 
resonant perturbation shifts all bifurcation points 
on both coexisting attractors in an opposite direc¬ 
tion so that one attractor is stabilized [Fig. 1(b), 
solution-1] while the second one is destabilized 
[Fig. 1(c), solution-2]. This means that by chang¬ 
ing the perturbation amplitude we may also con¬ 
trol at will the overlapping of different dynamical 
regimes associated with these new attractors. The 
third possibility to control overlapping is to change 
of the perturbation phase ip. 

Experimental and numerical examples of 
switching between new coexisting attractors when 
the system was initially in a 2 T state are shown 
in Figs. 2(a) and 2(b), respectively. For these ex¬ 
perimental conditions only one attractor exists in 
the phase space. After switching on the resonant 
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Fig. 2. Experimental (a) and numerical (b) stroboscopic 
CO 2 laser responses versus time (the laser intensity is sam¬ 
pled with the modulation period T). The unperturbed 
system (without the resonant perturbation) is initially in a 
2 T state, (a) The experimental parameters: Va = 3.35 V 
(p = 1.86), V p = 3.64 V, tp = 30°. (b) The parameters used 
in the computer simulation: /ro = 50 kHz, 7ro = 200 kHz, 
f d = 100 kHz, p = 2.1, <p “ 30°, fc d = 1.25 x 10" 4 , 
kp = 2.4 x 10~ s (fcp/fcd = 0.192), k s = 2 x 10 ~ 3 (k s /k d = 16); 
the point in time of turning on the resonant perturbation 
is shown by fl, the action of the pulse loss perturbation is 
shown by f2. 
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perturbation at t = 500 (units of the driving pe¬ 
riod T), the primary attractor splits into two ones. 
The system remains in the same 2 T state but with 
a slightly suppressed amplitude. Simultaneously, 
there appears a new attractor which is in a 4 T 
state. Acting by a single pulse perturbation (at 
t = 1000) the system was switched to the 4T state. 
The second pulse (at t = 1500) returned the sys¬ 
tem to the suppressed 2 T state. Both switchings 
were performed practically without transients. It 
is worth noting that the duration of the transient 
after the action of the pulse perturbation strongly 
depends on its amplitude and can be reduced for 
large enough amplitudes to one or two periods of 
the driving modulation. Similar results were ob¬ 
tained in numerical simulations which are in a good 
agreement with experimental ones [Fig. 2(b)]. 

A more interesting case is shown in Fig. 3 when 
the system was initially in a 4T state [Fig. 3(a)]. 
No other attractors were observed for these 
experimental conditions. Depending on its phase, 
resonant perturbations with the same amplitude 
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Fig. 3. Experimental stroboscopic CO 2 laser responses ver¬ 
sus time (the laser intensity is sampled with the modulation 
period T). The unperturbed system is initially in a 4T state. 
The arrow shows the action of the pulse perturbation. The 
experimental parameters: Vd = 4.11 (p = 2.28), V p = 3.52 V, 
(b) (p = 65°, (c) ip = 215°. 
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Fig. 4. Numerical stroboscopic CO 2 laser responses ver¬ 
sus time (the laser intensity is sampled with the modula¬ 
tion period T). The unperturbed system is initially in a 
4 T state. The symbol fl shows switching on the resonant 
perturbation, the symbol f2 shows the action of the pulse 
loss perturbation. The parameters used in the computer 
simulation: /ro = 50 kHz, 7 ro = 200 kHz, /a = 100 kHz, 
<p ^ 30°, k d = 1.55 x HT 4 (p = 2.6271), k p = 2.5 x 1CT 5 
(kp/kd = 0.1613), lj = 8x 10 -3 . 
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Fig. 5. Experimental stroboscopic CO 2 laser responses ver¬ 
sus time (the laser intensity is sampled with the modu¬ 
lation period T). The unperturbed system is initially in 
chaos. The experimental parameters: Vd = 4.84 (p = 2.69); 
(a) Vp = 0 V, (b) and (c) V p = 3.64 V, <p 9* 205°. 

may suppress periodic orbits [Fig. 3(b)] or induce 
chaos [Fig. 3(c)]. At the same time there appears 
a new attractor. The final state after switching 
also strongly depends on the phase of the resonant 
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perturbation. One can see that the system can be 
very quickly sent from a stable 2 T orbit to a stable 
8 T orbit [Fig. 3(b)], or from chaos to a stable 2 T 
orbit [Fig. 3(c)] and vice versa. It should be noted 
that the latter case [Fig. 3(c)] is very interesting 
from the standpoint of a fast change of the state 
in the system between “chaos” and “order”. The 
case represented in Fig. 3(c) was also simulated 
numerically and is shown in Fig. 4. It is seen to 
be in rather good agreement with the experimental 
results. 

Now let us consider the effect of a resonant 
perturbation in the bistability domain inherent in 
the system. Figure 5(a) demonstrates a coexistence 
of period-1 and period-3 attractors. The lat¬ 
ter appeared in the laser response after the ac¬ 
tion of the pulse perturbation. After switching 
on the periodic perturbation with the frequency 
/d/2, chaos of the period-! is converted to a 
2 T stable orbit, whereas the 3 T stable orbit is 
transformed to the 12T stable orbit [Fig. 5(b)]. 
This means that the resonant perturbation exerts 
a different effect on both coexisting attractors so 
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Fig. 6. Numerical stroboscopic'CO 2 laser responses versus 
time (the laser intensity is sampled with the modulation pe¬ 
riod T). The parameters used in the computer simulation: 
/ro = 50 kHz, 7ro = 200 kHz, /<j = 100 kHz, (p = 30°, 
k d = 2.15 x 10“ 4 (p = 3.6441), (b) k p = 2.5 x 10“ 5 ( k p /k d = 
0.11628), h 6 = 7 x 10“ 3 . 


that the period-1 attractor is stabilized while the 
period-3 attractor is destabilized. A different effect 
of the resonant perturbation makes it possible 
to control the overlapping of diverse states as¬ 
sociated with both coexisting attractors, and 
correspondingly to perform switchings between 
them. Similar results are obtained in the computer 
simulation shown in Fig. 6. 


4. Conclusions 

To conclude, we have presented experimental and 
numerical results on the dynamics of a CO 2 laser 
with modulated losses which demonstrate that 
resonant perturbations at subharmonic frequency 
not only radically change the dynamics of the sys¬ 
tem but also globally change the phase space by 
splitting the primary attractor into two new ones. 
We have demonstrated the possibility to control at 
will the overlapping of different dynamical regimes 
associated with these new coexisting attractors 
by changing the amplitude and the phase of the 
resonant perturbations as well as the bifurcation pa¬ 
rameter. We have also presented several examples 
showing that highly controllable switchings between 
these new attractors can easily be performed using 
the technique of short-lived large-amplitude pulse 
perturbations. Such a type of control of nonlinear 
systems could be of interest in the context of con¬ 
trol of systems with a fast response because it does 
not require any feedback system. The splitting of 
attractors and the switching between new attractors 
might be observable in a number of nonlinear sys¬ 
tems, both nonautonomous and autonomous, where 
the method of nonfeedback control of chaos has 
been implemented experimentally. 
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We demonstrate numerically and experimentally that a slow modulation of cavity detuning in a 
loss-modulated CO 2 laser can stabilize unstable periodic orbits even when the system remains 
in a particular dynamical regime for adiabatic changes of the detuning. When the parameter 
changes faster than the transient response of deformation of the original periodic attractor, the 
system can evolve toward an unstable periodic orbit. 


1. Introduction 

Recent trends in nonlinear dynamics are directed 
toward the development of effective methods of con¬ 
trolling chaos. A notable advance has been made by 
Ott et al. [1990] who proposed a feedback control 
technique based on stabilization of unstable peri¬ 
odic orbits embedded in a chaotic attractor. How¬ 
ever, in some kinds of systems (particularly, in bio¬ 
logical or chemical systems) a feedback loop is very 
difficult to realize. Lima and Pettini [1990] have 
shown that stabilization is possible without feed¬ 
back. They suggested the creation of stable peri¬ 
odic orbits from a chaotic system using a weak res¬ 
onant parametric perturbation. Both feedback and 
nonfeedback methods have been applied to many 
nonlinear systems including lasers (see e.g. [Roy 
et al., 1992; Bielawski et al., 1993; Meucci et al., 
1994]). 


‘Visiting from Stepanov Institute of Physics, Minsk, Belarus 


A simple nonfeedback and nonresonant con¬ 
trol method which does not require prior knowl¬ 
edge of the system behavior has recently been 
proposed by Yilaseca et al. [1996]. They have 
shown numerically that accurate stabilization of 
an unstable steady state in an autonomous sys¬ 
tem can be achieved by large-amplitude slow 
periodical modulation of a control parameter. 
Earlier, Liu and Rios Leite [1994] demonstrated 
numerically that stabilization of unstable periodic 
orbits in the Lorenz system can be achieved by 
coupling periodic modulation to a control param¬ 
eter. More recently, the main principles of the 
idea of Vilaseca et al. have been applied to a 
nonautonomous system and successfully realized 
in experiments with a loss-modulated CO 2 laser 
[Pisarchik et al., 1997]. The physical mecha¬ 
nism underlying the stabilization effect is tracking 
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unstable periodic orbits due to the delay of bifurca¬ 
tions when the control parameter is varied so that 
the system passes back and forth through an in¬ 
stability point [Mandel & Erneux, 1984; Kapral 
& Mandel, 1985]. A drawback of the method of 
Vilaseca et al. [1996] is the requirement for the am¬ 
plitude of the control modulation to be large (the 
larger the amplitude, the better the stabilization). 
As a consequence, the system response carries a 
large-amplitude envelope at the control frequency. 

In this work we present a method of stabilizing 
unstable orbits in a nonautonomous system which 
does not require very large amplitude modulation 
of a control parameter, as in previous studies with 
autonomous systems (see e.g. [Liu & Rios Leite, 
1994]). This new technique is similar to that pro¬ 
posed by Vilaseca et al. [1996] but the physical prin¬ 
ciple underlying the stabilization is quite different. 
We show numerically and experimentally that sta¬ 
bilization of unstable periodic orbits can be per¬ 
formed by periodic modulation of a system param¬ 
eter within some periodic or chaotic domain. In 
other words, to achieve the stabilizing effect it is 
not necessary for the system to pass back and forth 
through a bifurcation point. 

2. Model and Experimental 

Arrangements 

2.1. Model 

For numerical simulations of the operation of a CO 2 
laser with modulated losses we use a model based 
on the standard four-level scheme [Ciofini et al., 
1993]. Our modified model consists of eight dif¬ 
ferential equations for quasi-equilibrium and quasi¬ 
nonequilibrium populations of the upper and lower 
vibrational lasing levels, the global populations of 
the two manifolds of rotational levels, the popula¬ 
tion of the first excited level of N 2 , and the equation 
for average radiation density inside the cavity. 

The active medium of the CO 2 laser before 
switching on the electric discharge is a mixture of 
the gases CO 2 , N 2 and He. For the sake of definite¬ 
ness, we shall consider a single-mode lasing within 
a vibrational-rotational transition of the 00° 1- 
10°0 channel of the CO 2 molecule. According to 
Kuntsevich and Churakov [1994], the CO 2 laser 
with modulated losses can be described by the 


following system of equations 
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Here, 

Nq, N\ and N 2 are the 

relative quasi- 


equilibrium populations of the vibrational 00°0, 
10°0 and 00° 1 levels of CO 2 ; Mo and Mi are the 
relative populations of the fundamental and first 
exited vibrational level of N 2 ; 77 1 and 712 are the 
relative quasi-nonequilibrium populations of the vi¬ 
brational 10°0 and 00°1 levels of CO 2 1 ; n{ and ni, 
are the relative populations of lower and upper laser 
rotational sublevels 2 ; n e is the free electron density 
in active medium; W 2 \ and W\q are the effective 
rates of collisional relaxation in 00° 1-10°0 and 10°0- 
00°0 channels; Vr is the rotational relaxation rate; 
V\ and V 2 are the vibrational relaxation rates that 
describe the relaxation of “instantaneous” popula¬ 
tions 7i\ and n 2 to their quasi-equilibrium values, 
N\ and N2; Wnc is the exchange rate of the vi¬ 
brational excitation from N 2 to CO 2 ; /?i, @2 an d 
/?3 are the pumping rates of N 2 and the lower and 
upper levels of CO 2 in electric discharge; Nq and 


’The introduction of the populations Ni, N 2 , ni, 712 effectively allows us to take into account several important processes of 
intra- and intermode exchange in the active medium which are not included in Eqs. (1)—(3). 

‘For simplicity the rotational quantum number j is considered to be the same for both levels. 
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N n are the volume density of CO 2 and N 2 ; ko and 
A A: are the constant and additive losses; F{ and 
F% are the normalized Boltzmann functions deter¬ 
mining the part of molecules in the corresponding 
rotational sublevels in thermodynamic equilibrium: 
n\ = ^i n ii n 2 = -^ 2 n 2 i B( u ) an d K(t') are the 
Einstein coefficient and specific gain coefficient at 
the lasing frequency u and both have a Lorentzian 
lineshape; l m and l a are the lengths of the loss mod¬ 
ulator and active medium; c is the speed of light in 
the active medium; fi is the packing coefficient for 
the active medium in the cavity; u is the average 
radiation density; y — n 3 2 — n\ is the population 
inversion. 

The cavity loss coefficient is expressed as 
follows 

k = ko + -AA;[1 - cos(27 rf 0 t)] , (9) 

where ko is the constant losses without modula¬ 
tion, A k and fo = 1/T are the loss modulation 
(driving) depth and frequency, T is the period of the 
loss modulation. The expression for k is written in 
the form of (9) to satisfy the experimental situation, 
because the modulator of cavity losses used in our 
experiments always increases losses over their sta¬ 
tionary value ko, but does not decrease them. We 
define the cavity detuning as: 5 — (y — vo)h, where 
uo is the central frequency and 7 is the halfwidth of 
the Lorentzian gain lineshape. We introduce the 
control modulation of cavity detuning in a form 
similar to (9) 

S = So + ^A5[l - cos(27r/if)], ( 10 ) 

where 5o is the initial cavity detuning from the 
center of gain contour, i.e. without the control mod¬ 
ulation, A(5 and /1 are the amplitude and frequency 
of the control modulation. The modulation of S 
leads to the appropriate changes in the lasing fre¬ 
quency v. 

The numerical simulations are performed for 
the following set of parameters. The active medium 
consists of the gas mixture CC> 2 -N 2 -He=l-l -8 at 
the total pressure of 15 torr. The CO 2 laser oper¬ 
ates on a single mode at the 10P20 line. The cavity 
length is 2 m, the length of the active medium is 
1.8 m, ko = 6 x 10 -3 cm -1 , fo = 110 kHz. Other 
parameters (AA;, So, AS and / 1 ) are varied in nu¬ 
merical simulations. 


2.2. Experimental setup 

The experimental setup is similar to that described 
in previous work [Pisarchik et al., 1997]. The ex¬ 
periments have been performed on a single-mode 
CO 2 laser with modulated losses via an acousto¬ 
optic modulator. The driving electric signal, Vo = 
Ao sin(27r/ot), at frequency fo and amplitude To, 
is applied to the modulator providing the time- 
dependent cavity losses in accordance with Eq. (9). 
The control electric signal 

Vi=A\ + ^Ai [1 - cos(27r/i*)] (11) 

at frequency f\ and amplitudes of the constant A® 
and alternative A\ components, is utilized to tune 
the output mirror with the piezotranslator. This 
signal produces the appropriate changes in cavity 
detuning in accordance with Eq. (10). The fre¬ 
quency of relaxation oscillations at the center of the 
gain line estimated from averaged power spectra of 
the laser response to noise applied to the modu¬ 
lator is approximately 108 kHz. The frequency of 
the loss modulation, fo = 105 kHz, was chosen be¬ 
cause it matches one of the acoustic resonances of 
the modulator. This is close to the frequency of re¬ 
laxation oscillations at the center of the gain line 
but higher than the relaxation oscillation frequency 
in the region of significant detuning where the con¬ 
trol is established. Other parameters are the same 
as those used in the numerical simulations. 

3. Results and Analysis 

3.1. Diminution of periodicity 

Let us consider the effect of the modulation of cavity 
detuning when the CO 2 laser operates in a period- 
2 ( 2T ) regime. Figure 1 shows diagrams with 
the cavity detuning, So for the numerical diagram 
[Fig. 1(a)] and the voltage A® applied to the piezo¬ 
translator for the experimental diagram [Fig. 1(b)], 
used as a control parameter at relatively small val¬ 
ues of the driving amplitudes AA; and Ao- These 
diagrams are obtained without control modulation. 
Period-1 (T) and period-2 (2T) regimes are clearly 
seen on the bifurcation diagrams at certain detuning 
ranges. The windows of 2T behavior appear with 
detuning because the detuning moves the relaxation 
oscillation frequency of the laser into resonance with 
1/2T. 
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Fig. 1. (a) Numerical and (b) experimental bifurcation di¬ 

agrams with cavity detuning as a control parameter. The 
arrows indicate the range of detuning alternation within 
which the control will be applied, (a) A k = 1.4 x 10~ 5 cm -1 , 
A<5 = 0, (b) Ao — 7 V, Ai — 0. On the Y axis are plotted val¬ 
ues of the laser intensity sampled with the period T of the 
driving modulation. 


For stabilization of the unstable period-1 orbit 
that exists in the domain of the 2T regime, first, we 
select an initial point on the gain curve correspond¬ 
ing to the 2T regime (right peaks on the curves in 
Fig. 1). Then, following Eq. (10) we apply a sinu¬ 
soidal modulation to the cavity detuning so that the 
system remains inside the 2T domain. The ranges 
of the detuning variation are indicated in Fig. 1 by 
the arrows. 

Numerical (left) and experimental (right) stro¬ 
boscopic diagrams shown in Fig. 2 display the 
effect of the control modulation at different mod¬ 
ulation frequencies. One can see that with increas¬ 
ing control frequency /i, period 2 is progressively 


suppressed [Figs. 2(a)-2(c) and 2(e)-2(g)] and the 
unstable period 1 which is shown by dashed lines is 
stabilized as can be clearly seen in Figs. 2(d) and 
2(h). Comparing experimental and numerical re¬ 
sults one can see a good qualitative agreement. A 
difference in the shape of laser response between 
numerical and experimental curves is due to a dis¬ 
tinction in the contours of the numerical and exper¬ 
imental gain lines (Fig. 1) as well as because of a 
difference in the ranges within which the detuning 
is varied. 3 

3.2. Chaos suppression 

At high driving amplitude chaotic ranges appear 
in the bifurcation diagram at certain cavity detun¬ 
ings. To check the suitability of the method for the 
control of chaotic behavior of the laser, we choose 
a range of detuning modulation within one of the 
chaotic domains. Slow control modulation gives 
rise to the intervals of periodic behavior which al¬ 
ternate with chaotic ones at the control frequency 
[Figs. 3(a), 3(b) and 3(d), 3(e)] culminating with 
pulsed oscillations at the period 1 when fi further 
increases [Figs. 3(c)]. This effect is similar to that 
observed earlier in the experiments by Pisarchik 
et al. [1997] where the cavity detuning was larger 
and the system crossed T-2T bifurcation point. Al¬ 
though a complete stabilization of periodic-1 orbit 
is not achieved in experiment [Fig. 3(f)] because of 
noise, these results clearly demonstrate the validity 
of the present approach for inhibition of chaos. 

3.3. Discussion 

A possible physical mechanism underlying the ef¬ 
fect of dynamic stabilization can be associated with 
transient processes when the control parameter is 
changed. A detailed numerical and experimental 
study of transient processes occurring after fast 
change in the detuning within a period-2 domain 
allows us to reveal the following general features, 
(i) The duration of transient processes after for¬ 
ward and backward switching (with increasing and 
decreasing) the cavity detuning are different. When 
the cavity detuning increases, the transient pro¬ 
cesses are slower than when 5 decreases, (ii) When 
8 is switched from smaller to higher value, the sys¬ 
tem falls on an unstable period-1 orbit at any phase 


3 In the experiments the detuning does not pass through the maximal intensity at the 2T range. 
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Fig. 2. Time evolution of (a)-(d) numerical and (e)-(h) experimental sampled intensities, I, with period T showing destabi¬ 
lization of period 2 and stabilization of period-1 orbit with increasing frequency fi of control modulation within ranges shown 
by arrows in Fig. 1. (a), (e) /1 = 200 Hz; (b), (f) /1 = 500 Hz; (c), (g) /1 = 1 kHz; (d), (h) /1 = 2 kHz. Numerical parameters 
are So = 0.525, A 5 = 0.01. Experimental parameters are Aq = 7 V, Tfr = 5 V. Unstable period 1 is shown by dashed lines. 


of switching with respect to the loss modulation, 
and then evolves towards a stable period 2 with 
the leading Lyapunov exponent Ai = 7.7 x 10 2 s -1 . 
(iii) At the opposite switching (from higher to lower 
value), the system may also exhibit an unstable 
period-1 oscillation for a short time from which it 
later diverges with the leading Lyapunov exponent 
A 2 = 3.4 x 10 3 s _1 . The total duration of tran¬ 
sient processes after switching essentially depends 
on the switching phase, (iv) The dynamic stabi¬ 
lization takes place only when the period of the 
control modulation is shorter than 1/Ai. This corre¬ 
sponds to the control frequency f\ > 1 kHz. These 


simple estimations are in a good agreement with 
our numerical and experimental results presented 
in Secs. 3.1 and 3.2. 

All these features reveal a direct bearing of 
transient processes on the stabilization of unstable 
periodic orbits by periodic modulation of the detun¬ 
ing. The original periodic attractor is destabilized 
by the relatively fast change in the detuning. This 
suggests that the rate of change of the parameter 
requires a rate of change of the attractor that is 
greater than the negative Lyapunov exponent that 
governs the stability of the period-2 solution. Thus, 
the transition towards the unstable period-1 orbit 
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Fig. 3. (a)-(c) Numerical and (d)-(f) experimental stroboscopic (sampled with period T) diagrams showing inhibition of 

chaos with increasing control frequency f\. (a) f\ — 200 Hz, (b) fi = 1 kHz, (c) f\ = 6 kHz, (d) f\ = 200 Hz, (e) /i = 500 Hz, 
(f) fi = 3 kHz. Numerical parameters are So = 0.53, AS = 0.02. Experimental parameters are ylo = 10 V, A\ = 5 V. 


can be achieved when the relatively fast change in 
the parameter is faster than the transient response 
of deformation of the period-2 solution. 


4. Conclusions 

We have shown that dynamic stabilization of un¬ 
stable periodic orbits can be performed by a slow 
modulation of a control parameter even when the 
system does not pass through bifurcation points 
at adiabatic (quasi stationary) change of the pa¬ 
rameter. This method has been demonstrated nu¬ 
merically and experimentally in a loss-modulated 
CO 2 laser whose cavity detuning was periodically 
modulated. We suggest that the physical mecha¬ 
nism underlying the dynamic stabilization is con¬ 
nected with transient processes in the laser when 
the detuning is varied. Numerical simulations using 
an improved model of a loss-modulated CO 2 laser 
are in a good agreement with experimental results. 
Although the method has been implemented in a 
laser, we think that this is quite general and can be 


applied to many nonlinear systems whose initial 
state is characterized as chaotic or regular. 
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In this paper an examination is made of the performance characteristics of a scheme of chaos 
control using optimised impulsive delayed feedback. The scheme is applied in a model of 
external cavity laser diodes giving attention to the application of the feedback via modulation 
of the laser drive current taking into account practical constraints arising from technical delay 
and the frequency in application of the control signal. It is demonstrated that by use of a 
so-called e-ball condition chaos control is achievable in the fully developed coherence collapse 
regime without preliminary targeting dynamics to an unstable orbit. 


1. Introduction 

Interest in developing practical applications for 
chaotic dynamics has been stimulated by the 
development of techniques for controlling chaotic 
dynamics. The chaos control technique published 
by Ott et al. (OGY) [1990] has been successfully ap¬ 
plied to the control chaotic dynamics in a number of 
laser systems (see [Naumenko et al., 1998] for appro¬ 
priate references). In particular, the work includes 
attempts to control chaotic dynamics in external 
cavity laser diodes in the coherence collapse regime. 
Also mentioned is the use of occasionally propor¬ 
tional impulsive feedback (OPIF) [Gray et al., 1993; 
Liu & Ohtsubo, 1994], external microwave modula¬ 
tion [Ryan et al, 1994; Watanabe & Karaki, 1995; 
Liu et al, 1995] and continuous Pyragas-like elec¬ 
tronic feedback [Turovets et al, 1996, 1997]. These 
approaches have achieved rather moderate success 
due to specific features of semiconductor lasers such 


as high levels of intrinsic noise and the need for 
extremely high sampling frequencies. Interest in 
controlling semiconductor laser dynamics has been 
stimulated, in particular, by the possibilities for 
achieving secure communication systems which ex¬ 
ploit the properties of chaotic dynamical systems 
[Hayes et al., 1993]. An efficient algorithm for im¬ 
proving the locking rate between the receiver and 
transmitter of such a system has been reported pre¬ 
viously [Shore & Wright, 1994], The feasibility of 
optical injection-locking techniques for reciprocal 
synchronization of two distant chaotic laser diodes 
has been considered recently [Mirasso et al, 1996]. 
Another practical context where chaos control tech¬ 
niques may have an important role to play is in en¬ 
gineering immunity to coherence collapse caused by 
unintentional optical feedback w'hich may arise in 
the hybrid integration and packaging of commercial 
laser diodes where rather expensive optical isola¬ 
tors must, in general be used to avoid the relevant 
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optical feedback. Reviews of that effect and fur¬ 
ther references thereon are given in [van Tartwijk 
& Lenstra, 1995; Petermann, 1995]. 

The scheme of control of chaos proposed in 
the present study allows the convenient manipula¬ 
tion of the regimes of operation in such configu¬ 
rations. Alternative methods involving continuous 
optoelectronic feedback studied recently [Turovets 
et al., 1996, 1997] are not so flexible with respect 
to manipulating complex dynamics although they 
are possibly more advantageous in the context of 
engineering immunity to coherence collapse due to 
their modest requirements on the bandwidth of the 
electronic components in the feedback loop. 

With a view to such potential applications, this 
paper reports an investigation of the potential of 
discrete chaos control schemes based on the original 
OGY method [Ott et al., 1990] from its application 
to controlling the coherence collapse phenomenon 
in an external cavity laser diode. It is impor¬ 
tant to emphasize that the first attempts to control 
chaos in laser diode models were based on a over¬ 
simplified version of the OGY method like OPIF 
[Gray et al., 1993]. The need to readjust numer¬ 
ous control scheme parameters did not permit ef¬ 
fective optimization and lacked robustness to noise. 
For high operating frequencies the latter has had a 
rather discouraging impact on the perspectives of 

experimental implementation of such a scheme. In 
_ I 


contrast, in earlier work [Naumenko et al., 1998] 
attention has been focused on the optimization of 
parameters of control schemes with the aim to re¬ 
lax the requirements imposed by the basic features 
of laser diodes such as high levels of noise and fast 
dynamics. In particular, an examination has been 
made of several possible means of applying opti¬ 
mized discontinuous delayed feedback via different 
laser control parameters: the laser drive current, 
a laser field phase modulation in the external res¬ 
onator and modulation of cavity losses. In this pa¬ 
per we test this scheme using corrections applied to 
the laser bias current as the control without pre¬ 
liminary targeting to the unstable orbit -— which 
would require meeting the condition that two suc¬ 
cessive samples of laser intensity fall within a user- 
prescribed window: the so-called e-ball condition. 


2. Model 

It is assumed that the laser operates in a single 
longitudinal mode and is subject to weak optical 
feedback from an external mirror. A rate equa¬ 
tion treatment of this configuration has recently 
been developed on the basis of the Lang-Kobayashi 
model to take into account feedback effects in mod¬ 
ulated external cavity laser diodes [Langley et al., 
1995; Langley & Shore, 1993], The equations are as 
follows: 


dS(t') 

dt> 


dN(t) 

dt' 


Sit') 


+ ~* P + 2^ ^/S(fj^S(t' - Text) COS (0(f)) + F s (f) , 


m W) 


eF 


'sp 


- S(f)G n {N(f) - iV 0 }{ 1 + j + F n (t'), 


( 1 ) 

( 2 ) 


d ^p- = ± aGn {N(t')-N th } r- 

where 

9(f) = $(t') - $(<' - Text) + Wth^ext ■ (4) 

In the rate equations N(t') is the carrier den¬ 
sity, S(f) is the photon density, 3>(t') is the elec¬ 
tric field phase, 1(f) is the injection current and 
e is the electronic charge. Typical laser param¬ 
eters for a DFB laser are used, where V is the 
active region volume (1.5 x 10~ 16 m 3 ), r sp is the 


sin(9(i')) + F* (('), (3) 


carrier lifetime (2 ns), G n is the gain slope 
(2.125 x 1(T 12 m 3 s _1 ), N t h is the threshold car¬ 
rier density (9.9 x 10 23 m” 3 ), e' is the saturation 
parameter (3 x 10 -23 m 3 ), 7 is the spontaneous 
emission factor (1 x 10 -5 ), r p h is the photon life¬ 
time (2 ps), T is the confinement factor (0.4), ol is 
the linewidth enhancement factor (5.5), No is the 
transparency carrier density (4 x 10 23 m~ 3 ), w t h 
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is the operating frequency (A = 1.55 fi m), R 2 is 
the laser facet reflectivity (0.309), g is the laser to 
fibre coupling efficiency (0.4) and tl is the laser 
cavity round trip delay (9 ps). In the optical feed¬ 
back terms, i? ex t is the external reflector reflectiv¬ 
ity and T ext is the external cavity round trip delay 
(0.2 ns, which corresponds to an external cavity 
of approximately 2 cm for the fibre refractive in¬ 
dex 1.5). k ex t is the feedback coupling parameter, 
given by 


kext — 


1 VR^tv ■ 


The model can, in general, account for the ef¬ 
fects of Langevin noise terms, but the calculations 
described below consider noise-free deterministic 
dynamics. 

The introduction of new variables: u = SG n T sp , 
n = NG n Tp\\ , t = t'/ t s p , leads to the following nor¬ 
malized equations which are more convenient for 
the subsequent analysis: 


algorithm with a variable integration time step. 
Further technical details are provided in [Naumenko 
et al., 1998]. 


3. Impulsive Optoelectronic 
Feedback 

It is assumed that an impulsive electronic feedback 
signal is applied to the laser diode and that the am¬ 
plitude of every sequential impulse in the chain is 
proportional to the difference between the laser op¬ 
tical output powers taken at discrete moments of 
time. Such an arrangement implies detection and 
sampling of the output intensity synchronized with 
the pulse train generator. Mathematically, the im¬ 
pulse optoelectronic feedback terms in the dimen¬ 
sionless equations are taken in the form: 

K{t) =/3 + *i) 


u = vu(g(n, u)F — 1) + rn + 2 k^Juu r cos 6 , (5) 

h = Po + K(t) - n — ug(n, u), (6) 

4> = a(n — n t h)T — kJu T /u sin 0 , (7) 

6 = $ — $ r + wr. 

Here gin , u) = (n - n 0 )(l - £u), n 0 = N 0 G n T ph , 
£ = £'/G n T sp , W = WthTsp, r = I^T, V = T S p/Tp h , 

k = k ext T sp /T L , a = 0t'vT/2, P 0 = IG n T ph T sp /eV, 

T = Text/Tsp- 

The laser is taken to be biased by a DC injec¬ 
tion current Idc = 2 7 t h where 7 t h is the threshold 
current so as to ensure that the operation condi¬ 
tions are typical of the coherence collapse regime. 
Since the external cavity is short (2 cm), a typical 
route to chaos is a Hopf instability of a single exter¬ 
nal cavity mode steady state at R ext ~ 6.25 • 10~ 4 
or k « 2.5 followed by a Feigenbaum period dou¬ 
bling sequence. The first period doubling bifurca¬ 
tion takes place at A: « 4.26; full chaos is developed 
at k « 4.79. Referring to the classification scheme 
for regimes of operation of a laser diode with optical 
feedback (see e.g. [van Tartwijk & Lenstra, 1995] or 
[Petermann, 1995]), it is noted that due to the short 
external cavity the operating regimes II and III are 
mixed with the coherence collapse regime (IV) and 
the system enters regime IV directly from regime I. 

The terms K(t) in Eq. (6) describes an im¬ 
pulsive periodic control action applied to the driv¬ 
ing current. The rate equations (l)-(3) are solved 
numerically using a fourth-order Runge-Kutta 


- ll(lp + ti- Ti))f(t ~{(p + ti + T e l)) • 

( 8 ) 

The factor (3 is associated with amplification in the 
electronic feedback loop. The tunable phases ip 
specify a phase shift between a chain of applied cor¬ 
rection pulses and the reference signal which may 
be conveniently chosen as the laser output itself so 
as to measure phase with respect to a laser intensity 
maximum emitted at the discrete moments of time 
ti = J2k=i Fjfc. We will refer to the quantities /?, ip 
as the parameters of the control schemes which are 
to be optimized. Finally, Ti is assumed to be self- 
adjustable during the process of control time delay 
and equal to a time interval between successive laser 
intensity maxima. This means that the pulse gener¬ 
ator triggering time and the interval of sampling are 
synchronized to the detected maximum intensity. 

The term K(t ) is used to modify a selected 
control parameter of the laser diode slightly which, 
in the present case, is the driving current. Equa¬ 
tion (8) essentially assumes that the electronic com¬ 
ponents of feedback K(t) are fast enough to follow 
the laser output pulsations in real time, detect its 
maxima and trigger the pulse generator. In prac¬ 
tice, of course, there will always be some frequency 
cut-off due to the finite bandwidth and also techni¬ 
cal phase shift attributed in part to a finite speed 
of the electrical signal through the loop. This prac¬ 
tical aspect is taken into account via an additional 
time delay r e 1 in Eq. (8), which also includes a nat¬ 
ural time shift by a half-width of the pulse applied. 
It is supposed, nevertheless, that the bandwidth is 
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wide enough to permit microwave spectra of am¬ 
plitude fluctuations centered at the relaxation fre¬ 
quency (approximately 3 GHz for the parameters 
used here) to pass without significant distortion. 

The motivation for choosing feedback in the 
specified form is as follows. First, to affect stabi¬ 
lization the feedback should be negative, secondly, 
the feedback is designed so that it has no effect on 
the state which is supposed to be stabilized (e.g. an 
unstable periodic orbit with period T). Every kick 
of the control scheme has a functional form given 
by f(t). Its amplitude [see Eq. (8)] is proportional 
to the difference between the output powers that 
were emitted at instants of time ip + t* — r e \ and a 
previous time <p + U — r e i — T t . The functional form 
of the kicks was chosen to be super-Gaussian: 

f(t) = M exp(— at 2m ), a = (T(l/2 m)/mw) 2m , 

M — l/w , (9) 

where r(l/2?n) is the special T-function 
[Abramowitz &; Stegun, 1972]. 

The impulse function f(t) is normalized in such 
a way that its integral area is a unity. When its ef¬ 
fective width w tends to zero then it approaches a 
Dirac 5-function, whilst the exponent m regulates 
a specific shape of a kick. In particular, in the ex¬ 
treme case of m —> oo with a finite to, the pulses 
acquire a nearly square form. 

4. Results 

4.1. Targeting and control of 
unstable T-periodic orbits 

General strategies of control of chaos often imply a 
preliminary targeting to an unstable orbit which is 
to be stabilized. For the present system, which ex¬ 
hibits a period-doubling route to chaos, one possible 
way of achieving this is to pull a system back from 
a regime of developed chaos into a period doubling 
regime. In this way, the aim is to target and con¬ 
trol an unstable orbit and thence to track the sta¬ 
bilized system into the initial chaotic regime. Such 
an approach has been applied previously to low¬ 
dimensional dynamical systems such as modulated 
class B laser, and also to the external cavity laser 
diode configuration [Naumenko et al., 1998]. Here 
we present further optimization of such an approach 
with the aim to test the effectiveness of activat¬ 
ing the control scheme into the coherence collapse 
regime using an e-ball condition. 

First, the laser diode is assumed to be oper¬ 
ating with such a level of optical feedback that its 


dynamics is set just after the first period doubling 
bifurcation (k ss 4.6). The process of targeting con¬ 
sists of applying a single kick given by 

K^\t) = (3^ g f{t-p^ g ) ( 10 ) 

and applied at the appropriate phase <p t arg and 
with an appropriate amplitude /3targ> then after a 
short transient the system arrives in the vicinity 
of the unstable T-periodic orbit and finally, after 
staying a few periods on the unstable orbit the 
control scheme is activated. The difference signal 
(u(ip + ti) — u((p + U — Ti )) is constructed by sam¬ 
pling and storing the laser intensity at successive 
times and then a correction pulse is triggered and 
delivered to the driving current with a technical 
delay, r e [. 

In Figs. 1(a) and 1(b) the method is demon¬ 
strated. The first targeting impulse in this configu¬ 
ration is relatively large and is not shown here, but 
the moment of application of the targeting impulse 
is marked by a circle. It can be seen that the unsta¬ 
ble orbit is stabilized and the control signal tends to 
zero. Of course, the kick strength at the first stage 
depends strongly on how accurately the system has 
been targeted to the unstable T-periodic cycle. 

The effectiveness of the method is conveniently 
analyzed by computing domains of control and find¬ 
ing optimal conditions of control. In Figs. 2(a)-2(d) 
we present the domains of tracking the stabilized T 
orbit into the coherence collapse regime in the con¬ 
trol parameter space defined by the optical feedback 
strength k and the feedback technical delay r e i for 
different electronic feedback strengths (3 and fixed 
sampling phase tp. Computationally, this has been 
obtained by a successive stepwise changing of the 
optical feedback strength k or the technical delay r e i 
without readjusting the control scheme parameters 
at nodes of a grid in this parameter space and cal¬ 
culating the time, f re i, needed to relax once again to 
the T orbit with a prescribed accuracy (here taken 
to be 10 ~ 4 ) at every tracking step. Some small-scale 
details of the curves presented have been smoothed 
by spatial filtering. The external contour curve cor¬ 
responds to t re i = 300 and for practical purposes 
defines the external boundaries of the control do¬ 
main. Control is possible inside this domain and is 
not achievable outside of it. Along the inner curves 
in Fig. 2(a) t re i is equal to 100, 50, 35, 25 and 15 in T 
units. The curves in Figs. 2(b)-2(d) have a similar 
meaning, except that the curves with the smallest 
t re i are not shown. It can be seen that when an 
orbit has been stabilized one can vary the optical 
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Fig. 1. Time domain demonstration of the processes of targeting and subsequent activation of the control scheme to stabilize 
the unstable T-cycle using finite width kick perturbations of a driving current. The parameters of the control scheme are: 
normalized pulse-width w/T = 0.4; electronic feedback strength (j — 0.2; feedback phase <p/T = 0.58; form-factor of the 
superGaussian pulses m = 3. The parameters of the targeting pulse are ^targ/2T = 0.030825, /3targ = 0.1. Time is normalized 
to r sp , the normalized laser parameters are as follows: v — 1000, T = 0.4, r = 4 • 10~ 3 , no — 1.7, n t h = 4.2075, a = 2500, 
Ao = 8.2, ujt — 4.256, e = 7.059 • 10 -3 , r = 0.1, k = 4.6. (a) Laser intensity u versus time t. Having delivered the laser to the T- 
cycle by a targeting pulse [marked by the circle in Fig. 4(b)] the control scheme is activated [actually the first peak in Fig. 4(b)], 
(b) Controlling kicks signal versus time. Kicks measured in the same units as the pump current and tend to zero when control 
is achieved. 



Optical feedback strength 





Optical feedback strength 


Fig. 2. Domains of control in the control parameter space for the stabilizing technique using finite width pulse perturbation 
of the pump current with feedback-monitored by the difference of intensities. The solitary laser parameters are the same as 
in Fig. 1. Contour plot of the relaxation time to the T-periodic unstable orbit in the phase space: optical feedback strength 
k, versus technical delay r e i, the feedback phase ip = 0. The relaxation time t re i is measured from the moment when the laser 
has been prepared near the T-cycle up to the moment when the system relaxes to T-cycle with an accuracy of 10 -4 . The time 
is normalized to the period T. External curve corresponds to the f re i = 300 and gives, in practice, the external boundaries of 
the control domain. The inner curves are for t re i = 100, 50, 35, 25 and 15, respectively. The fixed electronic feedback strength 
/? = 0.2 (a), 0.3 (b), 0.4 (c), 0.5 (d). 




































































1796 A. V. Naumenko et al. 


feedback strength without modifying the parame¬ 
ters of the electronic feedback system and advance 
into the developed coherence collapse regime. It is 
noted that the relaxation times correspond to the 
reciprocal of the dominant Lyapunov exponent de¬ 
scribing the laser dynamics and hence the contour 
plots given in Fig. 2 provide an indication of the 
salient features of the dynamics. Moreover, the con¬ 
struction of these diagrams actually involves the op¬ 
timization of a particular control scheme, because 
the bottom of the landscape represents the control 
with fastest relaxation to an unstable orbit (and 
thus the state of the closed feedback loop with the 
largest negative Lyapunov exponent). In particular, 
it can be seen that there exists an optimal technical 
delay value centred around 0.25. This means that 
the unstable manifold of the stabilized orbit is ori¬ 
ented at these instants of time in such a way that 
a projection of the kick on it reaches the maximal 
value and so the result of control is the largest. It 
is worth emphasizing here that for optical feedback 
levels k ~ 5.9 the first large periodic window occurs 
in the bifurcation diagram of the system, so that 
the results presented in Fig. 2(b) show the possibil¬ 
ity of tracking the unstable T-periodic orbit through 
the entire first chaotic region. When the strength 
of electronic feedback is fixed at higher levels, it is 
possible to track the system even deeper into the co¬ 
herence collapse regime [Fig. 2(d)], however at such 
levels of electronic feedback the scheme is already 
unable to stabilize the T-orbit at k = 4.6. This 
implies that in order to track through the whole 
domain accessible by changing k it is also neces¬ 
sary to monitor and change some parameters of 
the feedback scheme (/? or p ). It is noted that, in 
accordance with its definition, the phase p of the 
applied correction kick is varied from 0 to 1 in units 
of T , i.e. the correction kick is applied every pe¬ 
riod. If, for example, the phase <p is more than 1 
(but less that 2), then the laser intensity is sampled 
and the correction pulse is applied to the driving 
current of every other period. Such a scheme, al¬ 
though still being formally described by Eq. (8), has 
different stability properties. In particular, the un¬ 
stable T- cycle in the Poincare section taken every 
other period has real positive Floquet multipliers 
(in contrast with the negative ones for the Poincare 
section taken every period) and the impulse control 
scheme described by Eq. (8) is not able to stabilize 
it. But the scheme is stable when the frequency 
of sampling is higher than the frequency of kicks 
application. So, the model presents the worse-case 



Fig. 3. Domains of control in the control parameter space 
the feedback strength /?, versus technical delay T e i, the feed¬ 
back phase tp = 0. The solitary laser parameters Eire the same 
as in Fig. 1. The normalized pulse-width w/T = 0.8 (b). 
The kicks applied every 2 T periods (solid lines), 3 T periods 
(dashed lines) and 4T periods (dotted lines). 

situation in the sense of the required speed of 
sampling of the laser intensity and applying the 
correction kicks. In fact, provided that the unsta¬ 
ble trajectory is subject to preliminary targeting 
and the optimal parameters of the control scheme 
are chosen, impulsive corrections can be applied 
less frequently (e.g. every second, third or fourth 
period). 

The effects of width and the application fre¬ 
quency of correcting pulses are demonstrated in 
Fig. 3 where the normalized pulse-width w/T = 0.8. 
Domains of control are shown in the control pa¬ 
rameter space of the feedback strength, /?, ver¬ 
sus the technical delay, r e i, at the fixed feedback 
phase p = 0. The kicks are applied every 2 T peri¬ 
ods (solid lines), 3 T periods (dashed lines) and 4 T 
periods (dotted lines). The sampling of the laser 
intensity is undertaken every period. It is seen that 
the domains of control basically retain their loca¬ 
tion and size, although their form is distorted as 
the frequency of application of the correction puls¬ 
ing is decreased. The stabilizing effect of increas¬ 
ing the correction pulse width is also clearly seen. 
Finally, the effect of alteration of the control do¬ 
mains with the period T (or recurrence with the 
period 2 T) with increasing the technical delay r e i 
may be also observed (see also [Naumenko et al., 
1998] where this effect is described and explained 
in more detail). It is expected that in the centres 
(optima) of the control domains the scheme also 
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offers the maximum robustness to noise. The results 
presented thus clearly show that such complex non¬ 
linear dynamic regimes as coherence collapse can be 
controlled with proper optimization. 

4 . 2 . Control of unstable orbits 
using an e-ball condition 

Attention is now given to the situation in which 
optimization of the control scheme, as identified 
in Sec. 4.1, has been made. In that situation it 


is possible to activate the control directly in the 
coherence collapse regime without preliminary tar¬ 
geting but by using a user prescribed window in 
which successive intensity maxima should fall — 
the so-called e-ball condition. Events of falling 
into the e-ball surrounding the unstable T-periodic 
orbit are guaranteed by the ergodicity of chaotic 
motion. This means that a chaotic trajectory must 
inevitably visit some neighborhood of any point em¬ 
bedded in the chaotic attractor. Of course, such 
events are stochastic by the definition and the time 
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Fig. 4. Time domain demonstration of the process of activation of the control scheme to stabilize the unstable T-cycle 
embedded in chaos area using the e-ball condition. The parameters of the control scheme are: the normalized pulse-width 
w/T = 0.8; the electronic feedback strength f3 = 0.4; the feedback phase <p/T = 0, the formfactor of the superGaussian 
pulses m = 3 and the technical delay r e i /T = 0.25. The parameters of laser diode are chosen in the coherence collapse zone, 
i.e. optical feedback strength k = 5 and the rest of parameters is as in Fig. 1. Time is normalized to r sp . (a) Maxima of laser 
intensity u versus time t. The control scheme is activated at time t = 0. After some chaotic transient the laser system finds 
itself in the vicinity of the unstable T-cycle satisfying to the e-ball condition, i.e. \u(ti) — u(tj_i)| < e = 0.05 followed by the 
short relaxation transient to the T periodic orbit stabilized, (b) Laser intensity u versus time t. (c) Control signal versus time. 
Measured in the same units as the pump current and tends to zero when control is achieved. 
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which elapses between activating the control scheme 
and achieving control depends strongly on the ini¬ 
tial conditions and should be treated statistically. 

In Figs. 4(a)-4(c) show how the scheme works 
in this case. The stroboscopic view (a), the laser 
intensity (b) and the control signal (c) are plotted 
as functions of time. Consideration is given here to 
the behavior of the laser after cessation of turn-on 
transient in the fully developed coherence collapse 
regime (with optical feedback strength k = 5.0). 
The control scheme was activated at time, t = 0, 
but the scheme was in fact in a waiting regime 
(/3 = 0) until the e-ball condition has been satisfied, 
i.e. (3 0 whenever \u((p+U)— T,)| < e and 

(3 = 0 otherwise. The first such event takes place at 
t = 5.5, but the very first kick of the control knocks 
out the system of the e-ball and the scheme is effec¬ 
tively switched off until the system next falls into 
the e-ball. After three such trials control is finally 
reached at t > 11. It is noted that the e-ball con¬ 
dition applied for control in Fig. 4 does not require 
a knowledge of the T-cycle at all. Thus, any occa¬ 
sional coincidence of two successive maxima turns 
on the control scheme. But it is still not guaranteed 
that this event takes place near a T-cycle, so some 
“false” alarms of the control scheme are possible as 
is shown in Fig. 4(c) (the first three kicks). We have 
also investigated a modification of this e-ball con¬ 
dition taken as \u(U) — u 0 | < e = 0.05, where uo is 
the maximum intensity of the T-cycle. It is found 
in this case that only one “false” alarm takes place 
at t = 1.5 and then at t > 3 control is successfully 
achieved. For this case it is found that on aver¬ 
age the stochastic transients in this case are shorter 
than those shown in Fig. 4. Finally, we observe 
that it is not required that uq be known a priori — 
rather it can be considered as a parameter of the 
control scheme. 


5. Conclusion 

In conclusion, a demonstration has been made of 
the opportunity for effecting control of the dynam¬ 
ics of external cavity laser diodes operating in the 
coherence collapse regime by the use of impulsive 
delayed feedback applied via the laser drive current. 
Account has been taken of practical constraints in 
the application of the feedback signal and control 
has been shown to be possible both with and with¬ 
out a preliminary targeting process. The effective¬ 
ness of the control scheme has been characterized 


via the computation of domains of control and by 
demonstrating the existence of optimal conditions 
of control. Such an optimized control configura¬ 
tion provides a significant advantage in practical 
application in comparison to previous chaos con¬ 
trol schemes which have been compromised by the 
noise properties of laser diodes. Moreover it has 
been shown to be possible to relax requirements in 
respect to the frequency of the application of the 
control corrections. From the viewpoint of practi¬ 
cal implementation this adds further to the attrac¬ 
tiveness of the present scheme by removing require¬ 
ments for very high frequency electronic circuitry 
for the application of the control. 
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An experimental investigation of the temporal dynamics of a diode laser subject to optical 
feedback from an external cavity containing a cell filled with cesium vapor has been performed. 
Peculiar dynamic regimes, such as self-pulsing, low-dimensional, and high-dimensional chaos, 
characterized by a new time scale, much longer than the time scales of all instabilities taking 
place in diode lasers under standard feedback conditions, have been identified. 


1. Introduction 

Diode lasers exposed to optical feedback from an ex¬ 
ternal cavity containing an atomic absorber are sys¬ 
tems of great interest in the fields of spectroscopy, 
atomic physics and quantum optics, as they pro¬ 
vide stable sources, with very narrow linewidth, on 
resonance with atomic transition frequencies. In re¬ 
cent years, these systems were widely investigated, 
mainly for the purpose of improving the spectral 
performance of the lasers and achieving an absolute 
frequency stabilization [Cuneo et al., 1994; Kitch- 
ing et al., 1994; Liu et al., 1994], On the contrary, 
there was no knowledge about the nonlinear dynam¬ 
ics and chaotic properties of these systems. In this 
paper, we provide a description of such properties 
and identify distinctive features of the dynamics re¬ 
lated to the presence of the absorber. 

We report on experimental investigations per¬ 
formed on an extended-cavity diode laser contain¬ 
ing a cell filled with cesium vapor, in a saturation 
spectroscopy configuration [Schmidt et al., 1994]. 


We focus on the temporal dynamics of the system, 
showing that the presence of the atomic absorber 
inside the feedback cavity produces characteristic 
dynamic regimes, which take place on a peculiar 
time scale, much longer than the time scales typi¬ 
cal of the dynamics of laser diodes under standard 
feedback conditions. Starting from a stable operat¬ 
ing regime characterized by self-locking of the laser 
frequency to the central Lamb dip of the saturated 
absorption spectrum of the Cs D 2 line, we have ob¬ 
served, with increasing detuning from the atomic 
resonance, a self-pulsing regime of the laser out¬ 
put in the form of a sequence of intensity dropouts 
and, eventually, a chaotic regime. Taking advan¬ 
tage of the fact that the dynamics of the system 
takes place on different time scales, spaced at least 
by three-four orders of magnitude, we have sepa¬ 
rated the slow dynamics induced by the absorber. 
Due to its long time scale, this dynamics is very 
easily accessible experimentally. The experimental 
time series have been analyzed with the time-delay 
method so as to reconstruct embedding portraits 
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and recurrence plots. The attractors in the phase 
space of the laser system have been characterized 
quantitatively by computing the Lyapunov expo¬ 
nent spectrum and the correlation dimension. From 
this analysis, it turns out that the signal observed 
for large detuning from the atomic resonance shows 
high dimensional chaos, characterized by two posi¬ 
tive Lyapunov exponents. This feature could either 
pertain to the hyperchaos [Rossler, 1979] or to the 
coexistence of different chaotic attractors. The dy¬ 
namics, either periodic or chaotic, observed out of 
the locking range, within a certain range of detuning 
from the atomic resonance, is characterized by the 
long time scale mentioned above. This time scale is 
absent in the dynamics of laser diodes in the stan¬ 
dard feedback set-up, and therefore, it identifies, 
among the operation regimes of the extended-cavity 
laser, those which are determined by the presence 
of the atomic absorber. 

Extended-cavity diode lasers with intracavity 
atomic absorber have never been studied previously 
from the point of view of nonlinear dynamics. In [di 
Teodoro et al, 1997], we discussed the spectral fea¬ 
tures of this system, pointing out that the presence 
of the absorber substantially modifies the mode pat¬ 
tern in the phase space of the extended-cavity laser. 
We also developed a rate equation model suitable 
for the interpretation of the spectral phenomena, 
which could not provide, however, a description of 
the global dynamics. In the absence of a compre¬ 
hensive theoretical model for the system considered, 
the analysis performed here makes it possible to 
infer general information on the dynamics in the 
phase space from the measured time dependence of 
the laser emission. 

This paper is organized as follows: In Sec. 2 
we describe the experimental apparatus and the 
measurements performed. In Sec. 3 we illustrate 
the observed dynamic evolution of the laser emis¬ 
sion and discuss the results of the time-series anal¬ 
ysis. Conclusions and final remarks are presented 
in Sec. 4. 

2. Experimental 

The experimental setup is shown schematically in 
Fig. 1. The observations were performed using 
a Spectra Diode Laser SDL-5400. The laser was 
optically coupled, through an antireflection coated 
collimating objective, to an external cavity termi¬ 
nated by a diffraction grating, placed 46 cm from 


the laser. A 4 cm long cell containing Cs vapor 
at the equilibrium density at room temperature 
was inserted within the external cavity. A frac¬ 
tion of the laser intensity was extracted from the 
external cavity by means of a beam splitter with 
less than 10% reflectivity and used for the detec¬ 
tion. To avoid unwanted feedback from optical 
surfaces in the measuring equipment, a magneto¬ 
optical isolator providing 40 dB of attenuation for 
counterpropagating light was placed on the path of 
the extracted beam. A scanning confocal Fabry- 
Perot spectrum analyzer, model 240 by Coherent, 
was used to detect the optical spectrum of the 
laser emission. Absolute optical frequency mea¬ 
surements were performed by comparing, with the 
spectrum analyzer, the laser frequency and the fre¬ 
quency of a reference diode laser. The reference 
laser was locked, through a standard technique of 
hybrid optical and electronic feedback, to the tran¬ 
sition F = 4 —» F' = 5 of the Cs D 2 line. The 
temporal evolution of the output laser intensity was 
detected, simultaneously with the spectral measure¬ 
ments, by means of an avalanche photodiode with 
1.2 GHz cutoff frequency. The photodiode signal 
was monitored with a Tektronix TDS540 digital os¬ 
cilloscope, with 500 MHz bandwidth. 



Fig. 1. Scheme of the experimental set-up. DL, diode laser; 
CO, collimating objective; Cs, cesium cell; BS, beam split¬ 
ters; G, diffraction grating; PZT, piezoelectric transducer; 
Iso, magneto-optical isolator; Ref, reference diode laser; FP, 
scanning confocal Fabry -Perot spectrum analyzer; APD, 
avalanche photodiode. 
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The diffraction grating terminating the exter¬ 
nal cavity was mounted on a piezoelectric trans¬ 
ducer, fixed to a mirror mounting. A continu¬ 
ous tuning of the laser wavelength within a range 
of tens of nanometers was achieved by combining 
the rotation of the grating and the variation of 
the cavity length through the piezoelectric trans¬ 
ducer. For laser frequencies far from the Doppler- 
broadened absorption line of the atoms, the feed¬ 
back power ratio was estimated, including losses at 
the optical surfaces within the external cavity, as 
about 30%. 

In the feedback scheme considered, the light 
field emitted from the output facet of the diode 
laser, tuned to the D 2 line, passes through the 
atomic vapor cell as a pump field and, after be¬ 
ing attenuated by the vapor, is retroreflected and 
used as a counterpropagating probe field, which is 
fed back into the laser cavity. At the center of any 
hyperfine transition of the D 2 absorption line, the 
pump beam bleaches a hole in the absorption spec¬ 
trum of the weaker probe beam. Therefore, the 
probe transmission, and hence the optical feedback 
intensity, is resonantly enhanced when the laser fre¬ 
quency matches one of the principal transition or 
crossover lines [Schmidt et al., 1994]. The satu¬ 
rated medium also modifies, due to its dispersive 
properties, the feedback phase. As an effect of the 
frequency-dependent changes produced by the in- 
tracavity absorber on the feedback parameters, we 
observed single-mode emission with frequency lock¬ 
ing to an atomic resonance, as well as multi-mode 
operation with each oscillation frequency locked to 
a different sub-Doppler feature of the atomic ab¬ 
sorption spectrum. This is a clear evidence of the 
influence of the atomic absorber on the laser oper¬ 
ation. An extensive analysis of the results of the 
spectral measurements performed on this system is 
reported in [di Teodoro et al., 1997], where a dis¬ 
cussion of the absorption and dispersion properties 
of the atomic vapor under saturation spectroscopy 
conditions is included. 


3. Temporal Dynamics and 
Time-Series Analysis 

Mono-mode self-locked operation was obtained 
within a locking range of ~160 MHz around the fre¬ 
quency of the hyperfine transition F = 4 —>• F' = 
5. The corresponding laser emission had a stable 
intensity, with a low noise level. Out of the lock¬ 


ing range, in the low-frequency direction, unsta¬ 
ble regimes in the emitted intensity were observed. 
The time dependence of the laser output was de¬ 
tected while varying, as a control parameter, the 
free spectral range (FSR) of the external cavity by 
means of the piezoelectric transducer. Here we in¬ 
dicate the relative difference of external-cavity FSR 
with respect to the case of self-locked operation 
as A v. 

Figure 2(a) illustrates a portion of the time 
series of the laser output intensity i(t ) recorded 
for Au =120 MHz. The laser output experiences 
dropouts which decrease the absolute intensity by 
about 30%. The dropout pattern appears to be 
periodic with a time scale of about 0.1 ms. The 
power spectrum of the signal, calculated by fast 
Fourier transform (FFT), and shown in Fig. 2(b), is 
characterized by a few, well separated dominant fea¬ 
tures. This is consistent with the picture of a pe¬ 
riodic signal. The observed period of the dropouts 
might be qualitatively associated with characteris¬ 
tic time scales of optical pumping processes within 
the absorber. For instance, the mean time for diffu¬ 
sion of Cs atoms at room temperature through the 
laser beam cross-section is of the order of tenths 
of millisecond [Schmidt et al ., 1994]. We should 
also notice that all documented intensity instabili¬ 
ties occurring in diode lasers subject to mere optical 
feedback from an external reflector take place on a 
much faster time scale. For example, in the regime 
of feedback intensity adopted in our experiment, 
extended-cavity diode lasers without absorber typ¬ 
ically exhibit dropouts of the output intensity, 
called low-frequency fluctuations (LFF), character¬ 
ized by a mean repetition rate of hundreds of MHz 
[Cerboneschi et al., 1993]. The intensity dropouts 
in Fig. 2(a) resemble, in shape, the LFF dropouts. 
However, in addiction to the difference in time 
scale, a clear distinction of the self-pulsing phe¬ 
nomenon detected in our system with respect to 
LFF is the periodicity. In fact, LFF appear as in¬ 
tensity fluctuations irregularly spaced in time and 
their phase-space dynamics is very complex and de¬ 
scribed in terms of chaotic itinerancy [Sano, 1994]. 
Sub-nanosecond dynamics, which is not detectable 
with our measuring equipment, is also involved in 
the low-frequency fluctuation regime [Fischer et al., 
1996]. 

The time series of the output laser intensity 
have been processed with the time-delay method 
[Packard et al., 1980]. Let i(t ) be the time de¬ 
pendent laser intensity, measured at equal sampling 
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Fig. 2. Laser intensity time-series (left-hand column) and corresponding power spectra (right-hand column): (a) and 
(b) Ai/ = 120 MHz, (c) and (d) Ai/ = 140 MHz, (e) and (f) Ai/ = 200 MHz. 


intervals S in the temporal window (0, r), where, for 
the signals recorded in our experiment, <5 = 2 /^.s and 
r = 32 ms. A portrait of the phase trajectory of the 
system is represented by the set of m-dimensional 
embedding vectors I„, defined as 

In = {i{n6), i(n5 + T ),..., i[nS + (m — 1 )T}}, 

n = 0,..., S, 

(1) 

where the lag parameter T is an integral multiple of 
the sampling interval 5 and S = [r — (m — l)T]/<5. In 
order to extract, from the time-delay portrait, both 
qualitative and quantitative information about the 
physics underlying the time series, a proper choice 
of m and T is needed. The algorithms applied in our 
analysis allowed us to test several values of m and 
verify the convergence of the calculated quantities 
for m > 3. The choice of parameter T was based 
upon the computation of the average displacement, 
according to the method outlined by Albano et al. 


[1988], Casdagli et al. [1991], and Rosenstein et al. 
[1994], The time-delay reconstruction was uti¬ 
lized to calculate recurrence plots, as introduced by 
Eckmann et al. [1987]. These plots are obtained, 
after choosing an embedding dimension m and a 
suitable radius r, as a square grid of S x S ele¬ 
ments, in which a dot is displayed at coordinate 
( k , n) whenever ||I^ —1^|| < r. The embedding por¬ 
trait corresponding to the time series in Fig. 2(a) is 
depicted in Fig. 3(a) and suggests that the dynamics 
of the laser intensity, in this regime, may be prop¬ 
erly described by a limit cycle. Consistently, the 
recurrence plot in Fig. 3(b) belongs to the so-called 
periodic typology [Eckmann et al., 1987], since its 
large-scale pattern is dominated by long lines par¬ 
allel to the main diagonal, which would be the only 
feature in the recurrence plot of a purely oscillatory 
signal, superimposed to an array of regular blocks 
emerging from stochastic fluctuations of the char¬ 
acteristic frequencies. 

A quantitative analysis of the time series has 
been performed by evaluating the largest Lyapunov 
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Fig. 3. Phase portraits for embedding dimension m = 3 (left-hand column) and recurrence plots for m = 10 (right-hand 
column) extracted from the experimental time-series: (a) and (b) A u = 120 MHz, (c) and (d) A u = 140 MHz, (e) and 
(f) A v = 200 MHz. The time lag used for the phase-portraits is T = 32, 24, and 16 ps in (a), (c) and (e), respectively. 


exponent, according to the method introduced by 
Rosenstein et al. [1993]. For the signal in Fig. 2(a), 
the largest Lyapunov exponent, calculated for dif¬ 
ferent embedding dimensions m, is nearly equal to 


zero, as expected for an asymptotically stable limit 
cycle [Parker & Chua, 1989]. 

In addition, the correlation dimension D corr 
[Grassberger & Procaccia, 1983] has been extracted 
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lnr 

Fig. 4. Logarithmic slope of the correlation sum C m (r) as 
a function of lnr, for several values of the embedding di¬ 
mension m: (a) Av = 120 MHz, (b) Ai / = 140 MHz, 
(c) Av = 200 MHz. The height of the plateaux indicates 
the value of the correlation dimension D c orr. 


from the experimental signals as 


Dr 


lim lim 

m —>oo r—>0 


d In Oui {v ) 
dlnr 


( 2 ) 


Here, the correlation sum C m (r) is given by 


Cm(r) 


2 

■S'(S'-i) 


0 ( r ~ 

{n,k} 



( 3 ) 


where 0 is the Heavyside step-function and {n, k} 
is a set of indices such that \k — n\ > t co /6, t co being 
a suitable cutoff time introduced to avoid artificial 
correlations arising from too closely spaced embed¬ 
ding vectors which may result from measurements 
taken nearly at the same time. In our analysis, t co 
was chosen as the inverse of the mean frequency 
of the power spectrum of the signal, obtained by 
FFT. The limit r -» 0 is, actually, unreliable, be¬ 
cause small values of r are blurred by noise and 
limitations on experimental accuracy. In practice, 
a plateau referred to as scaling region should ap¬ 
pear in the plot of d In C m (r) /dlnr versus lnr. In 
Fig. 4(a), such a plot is shown, for the time series 
in Fig. 2(a), for several values of m and the scaling 
region is clearly visible. The height of the plateau 
represents the estimated value of the correlation di¬ 
mension D corr . The result obtained, D con ~ 1, def¬ 
initely confirms that the dynamics is described by 
a limit cycle. 

In Fig. 2(c), part of the time series acquired for 
Av = 140 MHz is shown. In this case, the period 
fades in an irregular pattern of dropouts and the 
spectrum in Fig. 2(d) shows a broad-band struc¬ 
ture. Although broad-band spectra may be associ¬ 
ated with either stochastic or nonlinear determin¬ 
istic processes, the exponential decay of the spec¬ 
tral intensity at high-frequencies indicates chaotic 
dynamics [Brandstater & Swinney, 1987]. The dis¬ 
tortion of the phase portrait in Fig. 3(c) and the 
sharp change in the pattern of the recurrence plot 
in Fig. 3(d) support this interpretation. In par¬ 
ticular, in the recurrence plot, the long lines par¬ 
allel to the main diagonal appearing in Fig. 3(b) 
are now fragmented in shorter segments, the length 
of which has been proven to be inversely propor¬ 
tional to the largest Lyapunov exponent [Zbilut k, 
Webber, 1992]. A more striking evidence of chaotic 
behavior is provided by the largest Lyapunov expo¬ 
nent, which is now positive [see Fig. 5(a)] and by 
the value of D corr which is about 1.8 [see Fig. 4(b)]. 

The time series shown in Fig. 2(e) has been 
recorded for Av = 200 MHz. The signal shows 
no apparent regularity and the power spectrum 
in Fig. 2(f) has a broad structure which also in¬ 
cludes additional high-frequency components. The 
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Evolution Time ((is ) 

Fig. 5. Spectrum of Lyapunov exponents extracted from 
the time-series for m = 3: (a) Ais = 140 MHz, (b) A v = 
200 MHz. 

phase portrait in Fig. 3(e) is markedly different from 
the previous ones and the trajectory explodes in 
a blurred and folded geometry indicating that the 
embedding dimension m = 3 used for the portrait 
is too low to provide reliable topological informa¬ 
tion about the laser intensity dynamics. The recur¬ 
rence plot in Fig. 3(f) belongs to the homogeneous 
typology [Eckmann et al., 1987] and any small- 
scale texture is hardly visible. Again, a positive 
largest Lyapunov exponent is found [see Fig. 5(b)] 
and the calculated value of D COIr is about 6 [see 
Fig. 4(c)]. 

The characterization of the chaotic regimes has 
been completed with the computation of the entire 
spectrum of Lyapunov exponents. For this pur¬ 
pose, the algorithm introduced by Darbshyre and 
Broomhead [1996] has been implemented in a 


revised and noise-robust version described by 
Banbrook et al. [1997]. The procedure has been 
iterated in order to obtain the evolution of the 
Lyapunov exponents within a certain time inter¬ 
val and verify their convergence. The results are 
shown in Fig. 5, for an embedding dimension m = 3. 
Two positive exponents have been found for Av = 
200 MHz and only one for Av = 140 MHz. The 
calculation has been repeated for values of m up 
to 10, producing no substantial difference. Indeed, 
all additional Lyapunov exponents are negative in 
both cases. 


4. Summary and Conclusions 

The temporal evolution of the light emitted by a 
diode laser subject to feedback from an external 
cavity containing a cell of Cs vapor in a satura¬ 
tion spectroscopy configuration has been investi¬ 
gated experimentally. Frequency locking to the cen¬ 
tral Lamb dip of the saturated absorption spectrum 
of the D 2 line, accompanied by a stable laser emis¬ 
sion has been observed. With increasing detuning 
Av from the atomic resonance, within the range 
0 < Av < 140 MHz, where the laser operation 
appears to be strongly influenced by the presence 
of the atomic absorber, the output intensity dis¬ 
plays self-pulsing and low-dimensional chaos, which 
take place on a time scale much longer than the 
time scales typical of the dynamics of a diode laser 
coupled to an external cavity with no absorber. 
For larger values of Av, the laser emission shows 
high-dimensional chaos, characterized by two pos¬ 
itive Lyapunov exponents. Although the presence 
of two positive exponents is commonly associated 
with hyperchaos [Rossler, 1979], it might be ex¬ 
plained, here, in terms of coexistence of different 
attractors. The coexistence of attractors has been 
demonstrated theoretically also for diode lasers un¬ 
der usual feedback conditions [Masoller, 1994]. 

In the range of feedback strengths used in 
our experiment, high-dimensional chaos has been 
proven to underlie the dynamics of extended-cavity 
diode lasers with no absorber [Sano, 1994; Mirasso 
et al., 1997]. On the other hand, our analysis 
suggests that a low-dimensional attractor emerges 
when, for relatively small detuning, the laser oper¬ 
ation is dominated by the absorber. An interplay 
between different chaotic attractors might explain 
why no clear route to chaos has been observed. In 
diode lasers coupled to an external cavity with no 

































1808 F. di Teodoro et al. 


absorber, low-dimensional chaos, reached through 
definite routes, has been identified in several works, 
as, for instance, [M 0 rk et al., 1992; Ye et al., 1993]. 

The time-, frequency-, and optical frequency- 
domain analyses performed in this paper clearly de¬ 
tect the same aspects of the phenomena investigated 
and their mutual consistency has been regarded as 
a significant test. The variety of dynamic regimes 
identified in a range of laser frequencies close to the 
atomic resonance, along with the ease of the ex¬ 
perimental access to the slow absorber-induced dy¬ 
namics, makes this system an interesting subject for 
further experimental and theoretical investigations, 
which would be required for a full understanding of 
the nonlinear and chaotic properties described here. 
For instance, it would be needed to explain how the 
characteristic time scale induced by the absorber 
emerges from the destabilization of the regime of 
self-locking. 
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A very simple technique, which uses a succession of two different but constant levels in the 
pump parameter, is shown to be quite effective in controlling the turn-on transient of a class B 
laser. In particular, a drastic reduction in the delay time, as well as a lowering of the intensity 
overshoot under suitable conditions, has been observed both in numerical simulations and in 
experiments on a CO 2 laser. Such a straightforward, nonfeedback technique may have very 
promising applications in the realm of communications. 


1. Introduction 

Lasers, like many other optical systems, have at¬ 
tracted considerable attention for their nonlinear 
response to parameter variations and to perturba¬ 
tions. Their behavior, both regular and chaotic, 
was extensively investigated during the 1980’s 
[Abraham et al. , 1985; Bandy et al., 1988]. The the¬ 
oretical prediction that unstable periodic orbits in 
a developed chaotic regime could be stabilized with 
small perturbations [Ott et al., 1990] sparked in¬ 
vestigations in many different systems. Since then, 
much work has been done to control unstable pe¬ 
riodic orbits (e.g. [Hunt, 1991]) and in lasers sev¬ 
eral different schemes for the control of such orbits 
have been devised (e.g. [Roy et al., 1992; Bielawski 
et al., 1994]). More recently, interest has focused 
on reaching the desired periodic orbit in a minimum 
time, and many questions related to the targeting 
of particular states [Shinbrot et al., 1993; Boccaletti 
et al., 1997] have been examined. The application 
of optimal targeting to loss-modulated lasers has 
been very recently proposed in a few theoretical pa¬ 
pers [Kotomtseva et al., 1997; Turovets et al., 1997]. 
What is being stabilized in all these cases, however, 
is asymptotic behavior. Very recently, this field 


has expanded to address the characteristics of the 
transition between two states of a dynamical sys¬ 
tem and its control. A solution to the general 
problem of “customizing” a transient has been pi¬ 
oneered by Petrov and Showalter, who used a non¬ 
linear control method on a model for a nonlinear 
chemical reaction, with excellent results [Petrov & 
Showalter, 1996]. 

In this short communication, we present a sim¬ 
ple and effective way of controlling the transient 
evolution of a class B laser [Tredicce et al., 1985] be¬ 
tween two states: below threshold (which becomes 
unstable after the laser is switched on) and above 
threshold (final stable steady state). Although the 
control of laser turn-on may sound quite simple, 
class B 1 systems show nontrivial dynamical behav¬ 
ior due to the coupling of variables which respond 
at very different rates. Since a large part of all 
lasers sold for commercial applications (from high 
power to communications) belong to this class, a 
simple means of controlling their turn-on behav¬ 
ior could have very important and useful applica¬ 
tions. With the help of a simple model that cap¬ 
tures the essential dynamical properties of class B 
lasers, we discuss the basic physical features of these 
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1 The main elements of this class are solid state lasers (in particular Nd:YAG), semiconductor lasers, and CO 2 lasers. 
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systems when they are turned on by a sudden switch 
of the pump parameter. We then show the result 
of measurements obtained with a CO 2 laser and 
compare it to the predictions of a model that takes 
into account the specific molecular features of this 
laser. Finally we interpret physically the observa¬ 
tions and, in the conclusions, comment on possible 
applications. 

Extensive studies of the response of class A 
lasers (e.g. He-Ne, Ar + , and most lasers emitting 
in the visible part of the spectrum) have conclu¬ 
sively shown that their turn-on is controlled solely 
by the amount of spontaneous emission present at 
the switch-on time. The stochastic delay time with 
which these lasers respond by converting sponta¬ 
neous emission into a stimulated one closely mir¬ 
rors the statistics of the initial photons [Arecchi 
& Degiorgio, 1972; Arecchi et al., 1982] (below 
threshold). The response to noise at turn-on of 
class B lasers, which are modeled by two dynam¬ 
ical variables (instead of only one variable as in 
class A systems), has instead shown more complex 
features [Ciofini et al., 1992; Balestri et al., 1991; 
Balle et al., 1994; Grassi et al., 1994; Lippi et al., 
1997]. In this case, when the laser is turned on by 
abruptly increasing the pump value (which is what 
happens when a laser’s electrical switch is turned 
on), the system’s evolution is not just the result 
of a single variable, the intensity, that builds up 
from noise. Here the slow variable — the popula¬ 
tion inversion — is strongly coupled to the inten¬ 
sity and plays an important role in determining the 
system’s response, profoundly affecting its charac¬ 
teristics (cf. Sec. V in [Lippi et al., 1997]). 

The more complex response to noise in class 
B lasers is an indication of their intrinsic dynam¬ 
ical richness, which we exploit for the control of 
the turn-on transient. We are primarily interested 
in controlling the system’s switch-on on average, 
and so in the experiment we performed averages 
over many repetitions of the laser turn-on. In 
the theoretical discussion, we do not specifically 
include the influence of noise and instead explain 
our results using the deterministic elements of the 
model. Figure 1 shows the intensity of a CO 2 laser 
(lower trace) in response to a rapid variation of 
its pumping current (upper trace) from below to 
above threshold. Note the macroscopic delay time 
(~ 185 /is) between the abrupt current commuta¬ 
tion and the rise of the laser intensity out of noise. 
The overshoot in the intensity (past the final steady 



Fig. 1. Turn-on of a CO 2 laser by a sudden commutation 
of the pump parameter from below to above threshold. Up¬ 
per trace: Time behavior of the current flowing through the 
laser, which excites the molecules — initial current value: 
below threshold; final current value: above threshold. Lower 
trace: laser power. Note the delay time (« 185 gs) between 
the switch-up of the pump current and the actual laser’s 
turn-on, and the peak overshoot followed by strongly damped 
relaxation oscillations. 



tfiis) 

Fig. 2. Numerical integration of a five-dimensional model 
for a CO 2 laser, showing the intensity as a function of time 
after the application of an abrupt jump in the pump param¬ 
eter (see inset) at time t = 0. 

state value) is characteristic of these systems, as are 
the relaxation oscillations which follow. 

Figure 2 shows the result of integrating a 
complete five-dimensional model for a CO 2 laser 
[Meucci et al., 1992; Ciofini et al., 1993] for pa¬ 
rameters close to those of Fig. 1. This model, the 
best presently available, correctly reproduces most 
of the laser’s features: delay time, frequency and 
amplitude of the relaxation oscillations. Indeed, it 
was developed to account specifically for these fea¬ 
tures [Meucci et al., 1992] in the transient regime 
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and is therefore much more appropriate to the sim¬ 
ulation of the laser turn-on than equivalent reduced 
models (with two or three variables, e.g. [Oppo 
et al, 1989]). The improved agreement, also 
compared to earlier five equations models (e.g. 
[Arimondo et al, 1983]), is achieved using the re¬ 
laxation rates of the upper and the lower levels 
as fitting parameters. Although this procedure 
may seem somewhat arbitrary, it is justified by 
the fact that no simple dynamical model can take 
into account all of the laser’s features, especially 
in the transient regime (cf. [Meucci et al., 1992] 
for details). Our preference for this model comes 
from the fact that it predicts realistic delay times, 
while in the reduced models the delay is always 
largely underestimated (for reasonable parameter 
values). We remark that the time delay in the ex¬ 
periment (cf. Fig. 1) has an additional component 
of about 100 //s, which comes from the time re¬ 
quired to transfer the pump pulse to the lasing levels 
[Witteman, 1987] and which is never taken into ac¬ 
count in dynamical models. Therefore, the com¬ 
parison has to be drawn after subtraction of this 
additional time lag. In spite of its considerable 
advantages, one should not expect a quantitative 
agreement on the peak amplitudes from this model, 
either. Indeed, even in this case the amplitude 
of the peak overshoot is higher than in the real 
system. 

Although useful for the comparison between ex¬ 
periment and theory, the five-dimensional model is 
too complex for a simple analytical treatment, from 
which we would gain a physical understanding of 
the mechanisms involved in the transient. We will 
therefore discuss a reduced two-dimensional (rate 
equation) model [Siegman, 1986], which captures 
the main topological features of the phase space: 

- = (la) 

^ = - 7 ||[D(l + /)-P(()], (lb) 

where I and D are the e.m. field intensity and pop¬ 
ulation inversion variables, respectively, K and 7 || 
are their respective relaxation constants, and P(t) 
is the pump parameter, which is switched here at 
time t = 0 from below (P 0 ff) to above threshold 
(■Pon)- 


A numerical integration of these two equations 
is given in Fig. 3. In Fig. 3(a), the intensity re¬ 
sponse of the laser after the switch-on at time t = 0 
shows the effects of the reduced phase space in the 
increased number of relaxation oscillations present. 
The corresponding time evolution of the population 
inversion is shown in Fig. 3(b) where the horizontal 
dashed line at 1 represents the threshold value. The 
joint phase space portrait is given in Fig. 3(c). 

Since this model is purely deterministic, 
the below-threshold condition corresponds to the 
steady state (J Sj j = 0, D St i = P 0 ff)- Thus, the pop¬ 
ulation inversion D, as described by Eq. (lb), ini¬ 
tially undergoes an exponential growth and in the 
absence of laser intensity would asymptotically sat¬ 
urate at P on . The evolution of D closely follows 
this functional dependence as long as I is negligibly 
small: certainly whenever its growth rate [Eq. (la)] 
is negative, and in part even beyond that. As soon 
as threshold is reached, the sign of the prefactor 
of I in Eq. (la) changes [(1 — D) becomes nega¬ 
tive] and the intensity begins to feel the nonlinear 
(bilinear) coupling between itself and the popula¬ 
tion inversion. Thus, although the intrinsic growth 
rate of the intensity is much larger than that of the 
population inversion (K » 711 ), its actual growth 
depends on the evolution of D. In particular, we 
observe that the population inversion, D, must be 
larger than 1 [i.e. than the final steady state value 
(Isj — Pon, — 1 )] for the intensity to grow, 
and that the field intensity’s growth is proportional 
to the distance in phase space between the instan¬ 
taneous population inversion and its equilibrium 
value. 

Since the initial value of I is negligibly small 2 
(the spontaneous emission power is typically at least 
seven orders of magnitude smaller than the final 
laser power), there will be a certain time lag be¬ 
tween the moment when (1 — D ) becomes negative 
and when I assumes values that are non-negligible 
[compared to 1, cf. Eq. (lb)]. During this time D 
continues to grow [Fig. 3(b)], albeit at a slower rate 
[due to the term — y\\DI, Eq. (lb)]. The instanta¬ 
neous value of D, above D s j, causes the intensity I 
to overshoot and then to exhibit the characteristic 
damped, ringing behavior [Tredicce et al, 1985] well 
illustrated by the phase space plot [Fig. 3(c)] which 
clearly shows the various stages of the evolution of 
the transient. 


2 For the numerical integration we add a very small constant to I to allow the code to move away from the unstable fixed point. 
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D (norm, units) 

Fig. 3. Results of the numerical simulation of the laser rate equations under conditions comparable to those in Fig. 2. 
(a) Shows the laser intensity as a function of time after an abrupt commutation of the pump parameter at t = 0 (see inset of 
previous figure), (b) Shows the corresponding behavior of the population inversion, (c) Shows the phase space portrait. 


From the above discussion we extract the fol¬ 
lowing important points: 

1. There exists a deterministic component of the 
delay time at turn-on due to the slow response 
of the population inversion, which has to attain 
threshold before the intensity may grow. 

2. The field intensity cannot grow away from its ini¬ 
tial value until the population inversion is larger 
than 1 (in these normalized units). 

3. Due to the lack of a “macroscopic” intensity 
during the additional delay time, the popula¬ 
tion keeps growing further and further away from 
threshold under the steady action of the pump. 
(This “dead time” for the field is due to its rapid, 
but not instantaneous growth out of noise after 
the population is inverted.) 

If we want to influence the delay time with which 
the laser turns on, and, possibly also the height 


of the overshoot, we need to modify the evolution 
of the population. Since its relaxation constant yy 
cannot be changed, we must act on the pump pa¬ 
rameter directly. 

The simplest way of achieving our goal is that 
of adding to the pump pattern shown in the inset of 
Fig. 2 an intermediate control step (Fig. 4), whose 
height and duration can be chosen to suit our pur¬ 
poses. The addition of this intermediate step, of 
height P c (higher than the final level at P on ) and 
duration t c , serves a double purpose: (1) increasing 
the speed at which the population inversion grows 
— since its saturation value, P c , will be higher dur¬ 
ing this step — and (2) providing a better initial 
condition for the evolution of the field intensity (this 
point will be clarified in the physical discussion). 

Figure 5 compares the results obtained us¬ 
ing the five-dimensional CO 2 laser model [Meucci 
et al., 1992] for various values of the height, P c , 
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Fig. 4. Two-step control pattern for laser turn-on. The 
pump parameter is switched up from its initial value, P n n 
at time t = 0 , towards a control level of height P c , larger 
than Pan, and duration t c . P c and t c can be varied to op¬ 
timize the turn-on behavior. P on is the desired final pump 
level. 

of the control step. The effect of inserting this step 
is clear: The population inversion [Fig. 5(a)] in¬ 
creases more rapidly during the time interval during 
which the laser is pumped harder (between t — 0 
and t — t c — 20 /xs). As a result [Fig. 5(b)], 
the laser turns on sooner, i.e. the delay time at 
turn-on is reduced (mainly) by the faster growth 
of the population inversion from its initial value up 
to threshold. 3 A nontrival additional result, clearly 
visible in Fig. 5(b), is that there is an optimal height 
of the intermediate pump value for which we simul¬ 
taneously reduce the turn-on delay time and the 
peak intensity. 

The effect of the control step is summarized in 
Fig. 6, which shows the peak amplitude [Fig. 6(a)] 
and the delay time [Fig. 6(b)] as a function of con¬ 
trol step height. A strong deviation from a linear 
relationship (marked by the solid line in Fig. 6(b), 
obtained by interpolation between the leftmost two 
points) appears only after the minimum in the peak 
height has been passed [cf. Fig. 6(a)]. Near the min¬ 
imum the deviation from linearity is recognized to 
be approximately 8%. 

The experimental application of the pump pat¬ 
tern in Fig. 4 to a CO 2 laser (cf. [Grassi et al., 1994] 
for details about the setup), confirms the existence 
of both effects: The reduction in delay time and 
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Fig. 5. Numerical integration of a five-dimensional model 
for a CO 2 laser, (a) Temporal evolution of the population 
inversion for (a) A c = A on = 0.1, ( b) A c = 0 . 2 , (c) A c = 0.3, 

(d) A c = 0.4, (e) A c = 0.5, (/) A c = 0.6 and (. g ) A c = 0.7. 
The other parameters are t c = 20 ps, P G ff = 0.9, P on = 1.1, 
and the relaxation constants of the model k — 2 x 10 7 s _1 , 
71 = 8 x 10 4 s _1 , 72 = 1 x 10 4 s _1 , 7 n = 7 x 10 5 s _1 , z = 16. 
The horizontal dotted line (at 1) marks the lasing threshold. 
Notice that the duration of the control step corresponds to 
the time at which the population inversion curves (b-g) show 
a discontinuity given by the abrupt change in pump level, 
(b) Temporal evolution for the correspondent field intensity. 
The minimum in the intensity occurs between curves ( d ) and 

(e) , i.e. for a population inversion level at time t c very close 
to threshold (as predicted by the simplified model). 

the presence of a minimum overshoot height. In 
Fig. 7, we show the average values of peak height 
as a function of the average delay time for different 
amplitudes of P c (the values of the pump parame¬ 
ter are estimated under the assumption of a linear 
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3 The discontinuity appearing in Fig. 5(a) is not a numerical artefact, but reflects the sudden variation of the pump parameter 
onto the population inversion growth. We have extensively checked both the integration routine and the time step, in the 
framework of this and of another investigation [Lippi et al., 1997], to ensure that the sudden change in the control parameter 
does not affect the quality of the numerical prediction. 
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Fig. 6. (a) Height of the peak overshoot as a function of 

the control step; (b) time delay as a function of the control 
step; the straight line marks the extrapolated linear depen¬ 
dence of the delay, estimated from the two leftmost points. 
All parameters as in the previous figure (more A c levels are 
calculated here). 


relationship between excitation current and effec¬ 
tive pump). The averages are taken over 50 suc¬ 
cessive turn-ons. Three sequences, of 50 events 
each, are averaged together (for each value of 
the parameters) to provide the weighted averages 
with weighted standard deviations (cf. [Worthing & 
Geffner, 1943] for details), that we show in the fig¬ 
ure. The first measurement, far right point, marked 
as a, corresponds to P c = P on — i.e. turn-on with¬ 
out control (simple one-step commutation, as in 
Fig. 1). By gradually increasing the intermediate 
step height, we see that the average amplitude of the 
overshoot peak, plotted as a function of the average 



Fig. 7. Average peak intensity (weighted averages) as a 
function of average delay time (weighted averages) for dif¬ 
ferent values of the height of the control peak. Duration of 
the control t c = 20 ps. Estimated amplitudes of the control 
step: (a) P c « 1.25 (uncontrolled turn-on); (b) P c « 1.27; 
(c) P c « 1.28; (d) P c « 1.3; (e) P c « 1.34; (f) P c « 1.43; 
(g) P c « 1.61. P 0 fF « 0.9, P 0 n « 1.25. The intensity values 
are given without subtracting the detector’s offset, « 42 mV. 


delay time, decreases down to a minimum (« 2.5% 
reduction). If we continue increasing this ampli¬ 
tude beyond the optimal value, then the delay time 
is further reduced, but the overshoot amplitude in¬ 
creases. We first provide an approximate analytical 
explanation for the occurrence of this minimum and 
then interpret all results physically. 

We start from the equations of the simplified 
model already considered [Eqs. (1)]. For an analyt¬ 
ical treatment we use the approximation of [Balle 
et al ., 1991] where we neglect the intensity dur¬ 
ing the time period when it is very small. In 
this regime, 4 we can neglect the intensity term in 
Eq. (lb) and integrate the resulting expression for 
the population inversion as a function of time. Us¬ 
ing the equilibrium value of the population inversion 
D s> i as the initial condition, we obtain the value of 
the population at the end of the control step: 

D c = D(t c ) = D ofi e- Etc + P c ( 1 - e ~~ etc ), (2) 


where £ = y( 7 y /K), and t c is the time at which 
the control step ends. Using the latter expression 
[Eq. (2)] as a new initial condition, we integrate 
Eq. (lb) again, with pump value P on , to obtain the 


'As can be seen from Fig. 5 the minimum in the peak overshoot occurs when the pump is brought from P c to P on before the 
field intensity becomes macroscopic (a fact that will be clear later). Hence, the approximation we are presenting applies to 
the whole of the control step and beyond. 
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growth of the population inversion as a function of 
time after the control step. The approximate ex¬ 
pression that we obtain in this way is valid up to a 
time t*, which corresponds to the time for which 
the e.m. field intensity reaches a threshold value 
I*, which we can freely choose. Beyond this value, 
the intensity cannot be considered to be negligible 
anymore. We can write the value of the popula¬ 
tion inversion at time f as a function of all other 
parameters: 

D* = D(t*) = D oS ( 1 - st*) 

+ £[(1 + A on )t* + (A c — A on )£c] > (3) 

where 

P c = l + A c , (4) 

B 0 n = 1 + A on , (5) 

for convenience, and where we have expanded all 
exponentials to first order in st. 

After time t*, the full, coupled nonlinear system 
[Eqs. (1)] describes the dynamics. However, if we 
are only interested in determining the peak ampli¬ 
tude of the laser, we can once again solve the system 
of equations in an approximate form. In the region 
where the intensity is at its highest, the contribu¬ 
tion of the population inversion can, in comparison, 
be neglected [Balle et al., 1991]. Following the so¬ 
lution technique outlined in [Balle et al., 1991], we 
obtain for the peak intensity the expression: 

r, = r + + lo s D p~ io z D ‘ , (6) 

where I* and D* are the value of the intensity and 
of the population inversion at time t* , respectively, 
and I p and D p their corresponding values at the 
maximum of the intensity overshoot. If we now 
search for the existence of a minimum in the ampli¬ 
tude of this peak as a function of P c (or equivalently, 
A c ), we find that this occurs if and only if 

D c = l- A on s(t - t c ) , (7) 

which, because of the smallness of s, is very close 
to the laser threshold value (t is the time at which 
the population inversion reaches threshold during 
its growth). This is equivalent to saying that the 
minimum value of the overshoot is reached when 
the population inversion is brought very close to 
threshold by the control step. We note that the 


expression, Eq. (7), is independent of the “macro¬ 
scopic threshold” that we have chosen for the field 
intensity, I*. Looking back at the simulations in 
Fig. 5, we remark that the minimum in the inten¬ 
sity overshoot indeed occurs near the condition in 
Eq. (7). 

This result can also be understood intuitively. 
As already mentioned, the decrease in delay time 
is mainly due to the faster growth of the popula¬ 
tion inversion during the control step. However, if 
the control step is kept beyond the time for which 
the threshold is passed, the population inversion 
grows at a faster rate in the time interval where 
the field intensity grows. Hence, corresponding to 
the further reduction in delay time that is brought 
about by the longer control duration is a higher 
excess of population inversion when the field inten¬ 
sity reaches macroscopic levels, and hence a higher 
overshoot. This explains the increase of the peak 
amplitude when the optimal value for the mini¬ 
mum is passed. The existence of the minimum it¬ 
self can be understood in the following way. The 
amount of spontaneous emission present at all times 
determines the instantaneous value of the field in¬ 
tensity. Once the population inversion crosses the 
lasing threshold, the intensity starts growing from 
the instantaneous value that it has acquired at the 
time of crossing. A faster growth of the population 
inversion, below threshold, implies a faster growth 
of the spontaneous emission, and hence an initial 
condition for the intensity somewhat closer to the 
final steady state value I s j. Hence, provided that 
the control step does not continue above threshold 
— but that it actually stops just before reaching 
threshold [cf. Eq. (7)] — the faster growth of the 
population brings the intensity closer to the final 
operating condition, thereby reducing the time in¬ 
tervening between the crossing of threshold and the 
occurrence of the peak. The consequence of this is 
a (small) reduction of the overshoot height, accom¬ 
panied by a small contribution in the reduction of 
the delay time. 

Before concluding, we remark that the opti¬ 
mal value for the control to obtain the minimum 
overshoot can also be found by keeping the con¬ 
trol height, P c , constant and varying its duration 
t c . Indeed, if we approximate Eq. (2) to first order 
in et c an d equate it to the expression for the opti¬ 
mal value of D c , Eq. (7), under the further assump¬ 
tion that the laser is initially very close to threshold 
(DoS ~ l)j we obtain an approximate “area law” of 




1818 P. A. Porta et al. 


the form: 

A c t c =-—~ Aon(i-ic), (8) 

whose validity, however, is limited to short times 
(because of the approximation of the exponential 
time dependencies). This “area law” can be numer¬ 
ically verified to an acceptable degree even with the 
five-dimensional CO 2 laser model, but the equiva¬ 
lence between the time duration and the step height 
in the control exists only under very restrictive con¬ 
ditions. In general, a better degree of control is 
obtained by changing the step height, as we have 
found both numerically and experimentally. 

In conclusion, we have seen that a very sim¬ 
ple control technique, based on the introduction of 
an intermediate flat level in the pump parameter, 
is quite effective in reducing the class B laser’s de¬ 
lay time at turn-on and, under suitable conditions, 
even in lowering the height of the intensity over¬ 
shoot. The advantage of such a simple nonfeedback 
technique is that it could be implemented on very 
fast systems — e.g. semiconductor lasers — where, 
at present, any other technique of influencing tran¬ 
sients is unrealizable. In spite of the potential tech¬ 
nological difficulties that can be encountered in the 
implementation of even such a simple method, its 
applications could be rather far reaching. For in¬ 
stance, used to reduce the delay time at turn-on, 
it would permit faster modulation and higher bit 
transmission rates in VCSELs. 
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lier developments and discussions. P. A. Porta has 
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The synchronization of chaotic systems have received an increasing interest in the last few 
years. In an attempt to understand some of the possible mechanisms of synchronization of 
neurons in a noisy environment, the present study extends a control method which is based on 
the Kalman filter. This adaptive control mechanism is able to modulate the frequency and the 
clustering behavior of a network of neurons and can thus make the network switch dynamically 
to different rhythmic activity, leading to different coding possibility. 


1. Introduction 

Individual neurons often exhibit temporal chaotic 
behavior as observed in the characteristics of 
intracellular voltage measurements [Aihara &; 
Matsumoto, 1986; Mpitsos et al , 1988; Hayashi & 
Ishizuka, 1992; Abarbanel et al., 1996]. Individual 
neurons generate chaotic oscillations but groups 
of coupled neurons can display quasiperiodic 
synchronous rhythmic activity. Such dynam¬ 
ics are typical for many neurons in cortex and 
small neural systems like central pattern genera¬ 
tors that control the rhythmic motor behavior of 
animals [Eckhorn et al., 1988; Singer, 1993; Fujii 
et al., 1996]. Physiological studies have accu¬ 
mulated evidence indicating the existence of syn¬ 
chronous rhythmic activity in different areas of 
the brain of some animals [Gray et al., 1989; 
Engel et al, 1990; Bressler et al., 1993; Whit¬ 
tington et al., 1995; Neuenschwander & Singer, 
1996] and in some other small neural system 


[Dickinson et al., 1990; Meyrand et al., 1991; Wu 
et al., 1994]. For instance, based on in vitro ex¬ 
periments on the stomastogastric nervous system 
of a Crustacea, Meyrand et al. [1991] have reported 
that under an identified neuromodulatory stimu¬ 
lus, neurons operating independently as members 
of differents networks may be reconfigured into a 
new pattern-generating network that oscillates co¬ 
herently but differently from the original networks. 
Moreover, it is commonly assumed that every cog¬ 
nitive action is mediated by the coherent activity of 
neuron groups at multiple locations and it has been 
suggested [Eckhorn et al., 1988; Singer, 1993; Gray 
&; McCormick, 1996; Fujii et al., 1996; MacLeod 
&; Laurent, 1996] that this synchronous activity 
may have a role in solving the so-called binding 
problem. 

So some coordination or synchronization of 
these chaotic individual activities must be ar¬ 
ranged for a directed functional behavior that 
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have cooperative properties related to the func¬ 
tion performed. The cooperative behavior of these 
coupled neurons can be much richer and much 
more organized than the activity of the individ¬ 
ual neurons forming network. A fundamental ques¬ 
tion that remains open is how such assemblies are 
transiently self-synchronized for each specific task. 
In this view recent findings in the field of nonlin¬ 
ear dynamical systems are very promising, open¬ 
ing the possibility of controlling [Ott et al, 1990; 
Shinbrot et al, 1993; Chen & Dong, 1993] and 
synchonizing [Pecora & Carroll, 1990; Rulkov et al ., 
1995] chaotic dynamics. It has been shown that two 
identical chaotic systems can be synchronized by 
applying small temporal perturbations (parameter 
or variable perturbations) to one of them [Pyra- 
gas, 1992; Lai & Grebogi, 1993; Cazelles et al., 
1995], 

There has been some work on applying control 
techniques to networks of model neurons. Sepulchre 
and Babloyantz [1993] have used the Ott-Grebogi- 
Yorke control technique in a high dimensional net¬ 
work of coupled neurons and shown that higher 
dimensional control is possible. Lourengo and 
Babloyantz [1994] have used the Pyragas method to 
stabilize the unstable periodic orbits of the chaotic 
dynamics of network neurons of moderate size. A 
pulse control scheme that requires less knowledge 
of the system than an Ott-Grebogi-Yorke control 
technique is used by Carroll [1995a, 1995b] to con¬ 
trol and synchronize the dynamics of a network 
of four neurons described by a modified version of 
the FitzHugh-Nagumo equations. Carroll [1995b] 
has shown that for some coupling configuration the 
group of neurons are only partially synchronized 
even though all groups share a common drive. 
Giiemez and Matias [1996] have studied different 
coupling of small groups of simple chaotic neuron 
models which are synchronized by applying pro¬ 
portional perturbations to the dynamical variables 
of the system. Their method suppresses chaos 
in each neuron, yielding to a particular periodic 
or quasiperiodic dynamics. Abarbanel and his 
group [Abarbanel et al, 1996; Huerta et al., 1997] 
investigate the variation of the rhythmic activity 
produced by two coupled chaotic Hindmarsh-Rose 
neurons. They showed that sufficiently strong in¬ 
hibitory coupling between chaotic neurons organizes 
regular out-of-phase rhythmic behavior and that 
the change of the strength of the inhibitory coupling 
is responsible for the variations in the number of 
spikes in each burst. 


The purpose of the present work is to ex¬ 
tend a previous work [Cazelles et al., 1997] to 
the case of a network of model neurons. These 
previous works have pointed out that adaptive 
control even in noisy environment can provide a 
mechanism that permits a network of chaotic ele¬ 
ments to modify or stabilize its clustering behavior 
and can be necessary to maintain a network at a de¬ 
sired ordered state. These control mechanisms can 
also induce transition from high dimensional disor¬ 
dered motion to coherent synchronized pattern via 
feedback parameter perturbations. The adaptative 
control method [Cazelles et al, 1995] is based on 
the Extended Kalman Filter (EKF). This method 
progressively adapts the system parameters and/or 
the system variables to tune the dynamics towards 
a target behavior, taking into account uncertain¬ 
ties both in the system behavior and in the goal 
behavior. 

Depending on the parameters of the model net¬ 
work, different types of clustering behaviors from 
synchronized to completely desynchronized and dif¬ 
ferent types of dynamics from quasiperiodic to 
chaotic can be observed from a single network of 
neurons. In response to a certain external stimulus 
the network could modify its clustering behavior 
and the firing frequencies of the neurons. The 
modification of the rhythmic activity of the network 
is obtained by applying the EKF to force the neu¬ 
rons of the network to imitate a target dynamics, 
which could be the dynamics of a pacemaker cell 
which would want to impose its behavior on the 
network. Maintenance of an “optimal” function¬ 
ing of these neural systems depends on their ca¬ 
pacities of extracting pertinent information which 
is embedded in large noise. The key feature of this 
approach is to cope with this kind of noise to permit 
synchronization (for example) in spite of environ¬ 
mental changes and in the presence of large noise in 
the dynamics of the pacemaker. 

This article is organized as follows. In Sec. 2 the 
models used are presented and the EKF approach 
as a parametric adaptive control method is sum¬ 
marized. In the first part of Sec. 3 the dynamics 
of individual neurons are reported as the dynamics 
and the clustering behavior of the network. The 
second part of Sec. 3, after providing the conditions 
required to apply the control method, reports on 
how the degree of synchronization or the degree of 
chaos is modified as the parameters of the neurons 
are controlled. Finally, the last section is devoted 
to some concluding remarks. 
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2. Methods 
2.1. Models 

Actual neurons are quite complicated and for many 
years, simple approximate models of neurons have 
been used to understand their basic function. It 
has been shown in many cases that simple models 
can adequately reproduce some of the basic fea¬ 
tures behind the firing dynamics of neurons (see e.g. 
[Hoppensteadt, 1986]). The aim of the present work 
is to modify the synchronized dynamics of a network 
of neurons and also the firing frequency by adaptive 
control. So a simple model that reproduces the ba¬ 
sic features should be sufficient and in this present 
work, I used the Chialvo model [Chialvo & Apkar- 
ian, 1993; Chialvo, 1995] which has a neuron like 
dynamics and was found to agree qualitatively with 
experiments reported in previous works [Aihara & 
Matsumoto, 1986]. The model is a two-dimensional 
map written in the form: 

x(t + 1) = x(t ) 2 • exp (y(t) - x(t)) + k 
y(t + 1) = a ■ y(t) - b ■ x{t) + c 

where the x variable is related to an instantaneous 
membrane potential of the neuron and the y vari¬ 
able is equivalent to a recovery current. This model 
has four parameters, a determines the time con¬ 
stant of reactivation, b the activation dependence of 
the recovery process and c the maximum amplitude 
of the recovery current. The parameter k can 
be viewed either as a constant bias or as a time- 
dependent external stimulation. 

To simulate a network of neurons, each 
individual neuron is connected through a global 
coupling by the membrane potentials. One obtains 
globally coupled maps which are a dynamics of N 
individual cells evolving according to local map¬ 
pings and a “mean-field” interaction, the global in¬ 
formation influencing each individual cell. These 
maps are analogous to a “mean-field” version of 
coupled map lattices [Kaneko, 1989] in which the 
short-range coupling (nearest-neighbor diffusion) is 
replaced by a long-range coupling resulting from 
a feedback from the “mean-field” [Kaneko, 1992]. 
The model of the network of neurons is written as: 

j=N 

x(t) = (1 - e) • Xi (t ) + **(*) 

j=1 ( 2 ) 

Xi(t + 1) = x(f) 2 ■ exp (yi(t) — x(t)) + k 
y%{t + 1) = a ■ yi(t) - b ■ Xi(t ) + c 


where i is a time index, i a neuron index, N 
the number of neurons in the network, and s the 
coupling strength between neurons. 

2.2. The extended Kalman filter as an 
adaptive control method 

The basic idea of our control algorithm is to apply 
perturbations in the accessible system parameters 
to bring the system back to a desired state [Cazelles 
et al., 1995]. The adaptive control applies recur¬ 
sively to the system, an external force based on the 
difference between the actual output and a target 
behavior. This external force takes into account the 
correspondence between parameter values and be¬ 
havior types, the goal parameter value corresponds 
to the dynamics of the goal behavior. This con¬ 
trol method thus synchronizes the system dynamics 
with a target behavior. Various methods can be 
used for this task but the Kalman filter makes it 
possible to account for noise (uncertainties) both in 
the model and in the target dynamics. 

The evolution of the discrete nonlinear system 
(network of neurons in this study) is given by: 

X(f + 1) = /(X(<), p) + £(f) + External Force (3) 

X(t) is the n-vector of state variables, p is the j- 
vector of model parameters subject to control, / is 
a nonlinear map defined in (2) and £(f) the system 
noise (Gaussian noise with zero mean and variance 
matrix Q) The external force is computed by the 
Kalman control procedure. 

The target dynamics are produced by a 
reference model, a pacemaker cell that has the same 
dynamics as the individual neurons (1) but with 
constant parameter goal value p g . The underlying 
model of the target behavior is unknown for the sys¬ 
tem and only observations of the target behavior are 
available. An observation equation is defined: 

Y(t + l) = h(X(t),p g ) + V (t) (4) 

where Y (t) is the m- vector of observations, h the 
function which defines the relation between state 
variables (X) and observations (Y) and y(t) the 
observation noise (Gaussian noise with zero mean 
and variance matrix R) which can be interpreted 
as an error caused by a noisy environment. 

In (3) the accessible system parameter p 
operates as the control input driving the system to 
its target regime defined by noisy observations Y( t ).. 
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The change of the parameter vector y is defined 
according to: 

li(t + l) = /i(i) + f(t) (5) 

with £(f) being the parameter noise. 

The state equation (3) is increased by the 
expression (5) for the parameter dynamics and can 
be easily incorporated into a state-space framework 
as follows: 

Xa(t + l) = f(X a (t))+at) (6a) 

Y (t) = h(X a (t), fig) + Tf(t ) (6b) 


state-parameter noises; R the variance matrix as¬ 
sociated with the noise of the target behavior, $ 
and H are the jacobian matrices associated with 
the linearized form of the maps / and h ; the 
term K • [Y(f + 1) — h(X a (t + l|i))] represents 
the External Force (3) with K the Kalman gain 
matrix: 

K = P(f + 1|<) • H T (f + 1) 

• [H(f + 1) • P (t + l|i) ■ H T (f + 1) + Rf 1 

( 11 ) 

3. Results and Discussion 


where the increased state-parameter vector X a is 
defined as: Xj = [X J, p T ]. 

Originally, the Kalman filter is conceived for 
linear models and for nonlinear models a linearized 
system using the first order Taylor series expan¬ 
sion led to the “extended Kalman filter” (EKF) 
[Anderson & Moore, 1979]. EKF makes incremen¬ 
tal corrections of the state-parameter vector at time 
t to guide the system towards the target observa¬ 
tion vector of state-variables. For more biological 
realism, as in the work of Carroll [1995a, 1995b] only 
a pulse control has been performed. The corrections 
have been computed only when there is a spike in 
the target observations. The recursive form of the 
algorithm can be written as the following scheme: 

• Prediction of the states’ mean and variance: 

X a (t + l\t) = f(X a (t\t)), (7) 

P(t + l|f) = S(i+l|i)-P(i|i) 

•$(f + l|f) T + Q. (8) 

• Correction of the states’ mean and variance, only 
if there is a spike in the target dynamics: 

X a (t + l|f + 1) = X a (t + l\t) + K • [Y (t + 1) 

~h(X a (t + l\t)\, (9) 

P (t + l\t + 1) = [I - K • H(t)] • P (t + l|i) 

[I - K H(f)] T + K R K t . 

( 10 ) 

In the above equations, the notation (t\t) indicates 
an estimation at time t conditionally to all of the 
observations up to and including those available at 
t ; and superscript “t” denotes matrix transpose; P 
is the variance matrix of the state-parameter es¬ 
timation errors; Q is the variance matrix of the 


3.1. Neuron and network dynamics 

The individual neuron model used exhibits a rich 
dynamics, for example modifying the parame¬ 
ter b , one obtains either fixed point, periodic, 


(a) 



(b) 



(d) 




Times Index 


Fig. 1. Dynamics of the Chialvo model. Time evolution of 
the membrane potential with (a) b = 0.18, (b) b = 0.30, 
(c) b = 0.50 and (d) b = 0.60. Time series of the interspike 
intervals (e) obtained with different b values. 
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quasiperiodic or chaotic dynamics [Figs. l(a)-l(d)] 
[Chialvo, 1995; Guemez & Matias, 1996]. But more 
interesting from a neurophysiological point of view, 
when modifying the dynamics of the model the 
interspike interval is also modified; with b egal 
0.18 one has interspike interval chaotic time se¬ 
ries and with higher values one has fixed interspike 
intervals. But increasing the parameter b , the inter¬ 
spike intervals are increased, decreasing the firing 
frequency of the neuron [Fig. 1(e)]. It is particularly 
interesting because some neurobiologist conjectures 
that the firing frequency of a neuron acts as neu¬ 
ral coding [Fujii et al., 1996; MacLeod & Laurent, 
1996]. 

Chaotic globally coupled maps, as our net¬ 
work model, are known to have two conflicting 
trends: destruction of coherence due to the chaotic 
divergences of individual oscillator and a synchro¬ 
nization force through global averaging. In a glob¬ 
ally coupled chaotic system, the following clustering 
phases can appear [Kaneko, 1992]: (i) a coherent 
phase where all oscillators are synchronized; (ii) an 
ordered phase with few clusters in which each os¬ 
cillator is synchronized; (iii) a partially ordered 
phase with many coexisting clusters and (iv) a tur¬ 
bulent phase where all oscillators have their own 
dynamics. Figure 2 displays an overview of the 
clustering behavior and of the temporal dynamics 
of a network with 50 neurons (with 100 neurons 
quite similar results are observed). The clustering 
behavior diagram [Fig. 2(a)] is obtained by counting 
the numbers of different x values (membrane po¬ 
tential) for final attractor when 10000 points have 
been discharged. On one randomly chosen neuron, 
the temporal dynamics [Fig. 2(b)] is approximate 
by counting the number of different spikes for the 
last 5000 iterations. In function of the coupling 
strength, s, and of the parameter, 6, the behavior of 
the network varies from a phase of completely inco¬ 
herent chaotic behavior, through phases of partial 
synchronization to a global synchronization phase 
whereas the temporal dynamics can be chaotic, 
quasiperiodic or regular. 

Thus modifying the activation dependence of 
the recovery process, b, or the coupling strength, 
e, by adaptive control, one can induce different 
clustering behavior and different firing frequencies 
of the network thus modifying and arranging the 
coding possibilities of the network. 


3.2. Modification of the clustering 
behavior of the network 

Different types of clustering behaviors from 
synchronized to completely desynchronized and dif¬ 
ferent types of dynamics from quasiperiodic to 
chaotic can be obtained from a single network 
of neurons (Fig. 2). This opens the possibility of 
controlling the degree of chaos or the degree of 
synchronization of a network, by modifying either 
the parameters of each individual neuron or the 
coupling strength among neurons. In this work, 
I focus on modifying the activation dependence of 
the recovery process, b, by applying the Kalman 
control process described above. To use the EKF 
procedure, one must specify the target dynamics 
and matrices Q and R: (i) An external signal, the 
dynamics of a pacemaker cell which would want to 
impose its behavior on the network is used for ob¬ 
servations (the target behavior); (ii) the variance 
matrix Q of the state-parameter noises will take 
on values between 5% and 10% of the mean of the 
network dynamic; (iii) the variance matrix of the 
observation error R will take on the same values as 
Q, in the absence of noise; if observations (the tar¬ 
get behavior) are corrupted by noise, the variance 
of the additive noise will be assigned to R. 

At one time, I chose to synchronize an initially 
chaotic and turbulent network of 50 identical neu¬ 
rons (b = 0.18) with low coupling strength (s = 
0.02) to the dynamics of a pacemaker characterized 
by a quasiperiodic behavior ( b g = 0.25). Figure 3 
illustrates this example: Fig. 3(a) displays the pa¬ 
rameter adaptation; Fig. 3(b) shows the dynamics 
of the 50 neurons and Fig. 3(c) the decimal 
logarithm of the difference between the membrane 
potential of one randomly chosen neuron and the 
membrane potential of the pacemaker cell. The net¬ 
work and the pacemaker cell dynamics originally 
behaved quite differently. As soon as the control 
action is switched on at time index egal 100, the 
parameter value b starts converging towards the 
goal value b g [Fig. 3(a)] and the trajectories be¬ 
come fully synchronized [Fig. 3(b)], the difference 
between the two signals (network and pacemaker) 
decreased dramatically [Fig. 3(c)]. At time in¬ 
dex 400 the target behavior is modified, the new 
driven pacemaker cell have a quasiperiodic behavior 
with low firing frequency ( b g = 0.50). Once 
again the network is quickly resynchronized to the 
new rhythmic activity by the adaptive controller. 
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Fig. 2. Phase diagram for a network of 50 neurons. Part (a) shows the clustering behavior and part (b) the temporal 
dynamics of the network. The numbers of clusters, k, have been obtained by counting for the final attractor the numbers of 
different x values (membrane potential) when 10000 points have been discharged. Temporal dynamics of the network have 
been approximated by counting on one randomly chosen neuron the number of different spikes for the last 5000 iterations. All 
calculus were carried out by incrementing the parameters b and s by 0.001. 























Synchronization of a Network of Chaotic Neurons 1827 




0 200 400 600 800 

Time Index 
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Fig. 3. Synchronization of a chaotic turbulent network of identical neurons (N = 50 , b — 0.18 and e = 0.02) to the dynamics 
of a pacemaker cell characterized by a quasiperiodic behavior ( b g = 0.25 and b g — 0.50). The control is active at time index 
100 and at time index 400 the dynamics of the pacemaker cell is modified, (a) Parameter adaptation (blue fines: parameter of 
each neuron; red line: parameter of the pacemaker cell which defined the target dynamics), (b) Time series of the 50 neurons, 
(c) Time evolution of the log of the absolute value of the difference between the membrane potential ( X ) of one randomly 
chosen neuron and the membrane potential (X goa i) of the pacemaker cell. 


At each time, after 3 or 4 pacemaker spikes the pa¬ 
rameter convergence is complete and the dynamics 
of the network is fully synchronized to the rhythmic 
activity of the pacemaker cell. 

The robustness of the method used when noise 
is added to the target behavior have been tested. 
The noise intensity is characterized by the signal to 
noise ratio (SNR) which is the standard-deviation of 
the reference dynamics (noise-free target behavior) 
divided by the standard-deviation of noise. Figure 4 
displays the parameter adaptation and the be¬ 
havior of an initially chaotic and turbulent net¬ 
work (b = 0.18, s = 0.02) that is control to a 
target behavior with a quasiperiodic synchronized 


behavior ( b g = 0.25). For each neuron of the 
network, the dynamics of the pacemaker cell is 
spoiled by large additive noise (SNR = 5.00). As 
in the previous noiseless example, despite the large 
noise, the control successfully drives the parame¬ 
ter b towards the goal value b g and the network 
presents a synchronous rhythmic activity [Fig. 4(c)]. 
Figure 4(a) displays the parameter adaptation, 
Figs. 4(b) and 4(c) show the networks’ behavior 
before and after controlling. 

Similar results (not showed) are obtained with 
other values of the coupling strength (e) or the 
activation dependence parameter (6), for exam¬ 
ple with chaotic synchronized or with two-cluster 
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(a) 




Fig. 4. Synchronization of a chaotic turbulent network of identical neurons (N — 50, h = 0.18 and e = 0.02) to the dynamics 
of a pacemaker cell characterized by a quasiperiodic behavior ( b g = 0.25). For each neuron, the target behavior is spoiled by 
a large additive white noise characterized by SNR = 5.00 (7.0 dB). The control is active at time index 100. (a) Parameter 
adaptation (blue lines: parameter of each neuron; red line: parameter of the pacemaker cell), (b) Plot of the evolution of the 
first neuron of the network versus the 49 others before activation of the adaptive control, (c) Plot of the evolution of the first 
neuron of the network versus the 49 others after activation of the adaptive control (simulations were run over 400 iterations 
and the last 200 points are shown). 


periodic target behavior. The parameter fluctua¬ 
tions produced by the adaptive controller stabilize 
the network momentarily into a coherent state with 
a given firing frequency, more appropriate for a spe¬ 
cific information task. This behavior may explain 
short episodes of synchronization of the assembly 
of neurons at different firing frequency even when 
the external signal is embedded by large noise. At 
each time, synchronizations arise only from the ac¬ 
tion of an adaptive controller driven by an external 
input, the pacemaker cell. These synchronizations 
occur in sufficient short time scales that it may be 
reasonable to think that they could be explained by 
electrophysiological mechanisms. 


4. Conclusion 

Neuronal sensory systems of living organisms are 
capable of dealing with noisy signals. An impor¬ 
tant issue in neurosciences is the study of the abil¬ 
ity of neuronal systems to cope with uncertainty 
and the lack of precision which arises from different 
sources of perturbations. Temporal coherence or 
synchronous firing, postulated as a mechanism for 
functional organization of the assembly of neurons 
and as a potential neural code. 

In an attempt to understand some of the pos¬ 
sible mechanisms of synchronization of neurons in 
a noisy environment, this paper reports on an 
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adaptive control mechanism which is able to modu¬ 
late the frequency and the clustering behavior of an 
assembly of neurons. The adaptive control method 
used is analogous to a learning process: by com¬ 
paring the system dynamics with a target behavior 
and by correcting its dynamics, the system learns 
to adapt itself in so as to recover the target dy¬ 
namics, which could be the dynamics of a pace¬ 
maker cell. Adaptive control of coupled neurons 
may be an analogous mechanism to mediate the 
general process which is responsible for dynamical 
“linking” and “unlinking” of neurons into varying 
coherent groups. This controller can thus make the 
assembly of neurons switch dynamically to different 
rhythmic activity (e.g. from noncoherent to 
synchronized, from chaotic to periodic), leading to 
different coding possibility. 

The framework of control and synchronization 
of chaos seems promising for the treatment for 
the description of rapid changes in the dynamics 
of neuronal networks. If the assembly of neurons 
could be viewed as a network with spatiotempo- 
ral chaotic activity then rapid transition between 
states of different synchronized activity could be a 
framework that enables us to propose a mechanism 
for understanding some aspects of neural dynamics 
and recent works go along these lines [Lourengo & 
Babloyantz, 1994; Caroll, 1995a, 1995b; Giiemez 
& Matias, 1996; Abarbanel et al., 1996; Huerta 
et al., 1997]. 
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The analytical properties of the solution of a system of ODEs in the complex time plane 
influence its dynamical behavior on the real time axis. In particular, the extrema of the real 
time solution can be associated to the singularities of the complex solution falling close to the 
real time axis. Moreover for a twice differentiable stochastic process, the expected value of the 
number of extrema for unit time can be determined. These two results are used here as the 
starting point to introduce two new algorithms to test for time series nonlinearity. They do 
not require the phase space reconstruction protocol and seem to work well also for short data 
sets. 


1. Introduction 

The nonlinear analysis of time series is based on the 
theory of dynamical systems: It is assumed that 
the signal under examination is generated by a 
deterministic law. When experimental data are 
used, the analysis of time series becomes more dif¬ 
ficult for the presence of noise. However it is gener¬ 
ally assumed that the stochastic component in the 
signal is small and that it does not destroy com¬ 
pletely the original properties of the bare (without 
stochastic perturbation) dynamical system. Since 
a nonlinear system can exhibit deterministic chaos 
a first step, in the analysis of an irregular signal, 
is to check whether its values are characterized 
by nonlinear time correlations. Understanding the 
nature of the dynamical processes responsible for 
the observed data is one of the most challenging and 
important problems in nonlinear time series analy¬ 
sis and several methods have been proposed to this 
aim [Abarbanel et al. , 1993]. Here we suggest the 
use of two new algorithms to detect nonlinearity in a 


stationary time series (that is to check whether non¬ 
linear time correlations are present among the time 
series values) both based on the analysis of the 
extrema (local maxima and minima). Two theo¬ 
retical considerations justify these algorithms. The 
first is connected to the theoretical and numerical 
results obtained in the study of systems of ordinary 
differential equations (ODEs). It has been shown 
that the dynamical behavior of the real time solu¬ 
tion of an ODE is strongly connected to its analytic 
properties in the complex time plane, and in partic¬ 
ular to the distribution of the singularities nearest 
to the real axis [Ramani et al., 1989]. The second 
consideration arises from a general property of a 
stochastic process which states that given a mean- 
square differentiable stochastic process x(t) the 
expected number of its extrema for unit time is con¬ 
tained in the joint density function of x(t), x(t) and 
x(t) [Soong, 1973]. These theoretical and numerical 
results suggest that the sequence of the extrema of 
a time series contains dynamical information on the 
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process responsible for its generation. The method 
of surrogate data was used to test if a suitable statis¬ 
tic evaluated on a sequence of extrema could be 
accounted for by a Gaussian linear stochastic pro¬ 
cess [Theiler et al, 1992]. 

2. The Algorithms 
2.1. PSC algorithm 

The problem of establishing whether a system of 
ODEs is integrable or chaotic is generally diffi¬ 
cult. During the last two decades some impor¬ 
tant theoretical results have been obtained: Of par¬ 
ticular interest, a strong connection between the 
dynamical properties of the real time solution of an 
ODE and the analytical ones in the complex time 
plane, has been shown to exist [Ramani et al, 1989]. 
The solutions of the Lorenz equations in the com¬ 
plex time plane, for different dynamical regimes, 
were studied both analytically and numerically 
[Tabor & Weiss, 1981]. It was found that the dis¬ 
tribution of the singularities in the complex time 
plane is critical in determining the behavior of the 
real time solution. In particular, from the func¬ 
tion theory of complex variables, it follows that the 
main contributions to the real time solution come 
from the nearest singularities to the real time axis. 
This last fact is justified by the Cauchy’s integral 
formula for a function of complex variable /(z) an¬ 
alytic and single-valued in a simply connected do¬ 
main D: ip /(z)/(z— zo)dz = 2mf(zo) where T € D 
and zo is enclosed whithin T. For a better clarifi¬ 
cation we just consider the example of the Duffing 
equation x + ax + x 3 = bcos(t ) + c where a, b, c are 
real constants. Let us consider initially the very 
simple case a = b = c = 0: x + x 3 = 0; we search 
for a solution of this equation over the complex time 
plane. It is easy to prove that a particular solu¬ 
tion of this equation can be expressed by means 
of an elliptic function: x p (z) = cn(z,l/2), where 
z = t + it] x p (z) is a double periodic function of z. 
In general the function cn(z, q), q 2 < 1, is analytic 
over the complex time plane except for the set of 
points A = {z m>71 = 2mK + i(2n + 1)K' : m,n € Z} 
where polar singularities are present; K and K' are, 
respectively, the complete elliptic integral and its 
complementary (constants for a fixed q ). From the 
definition of the set A it follows that the nearest 
singularities corresponding to the real time axis are 
those for which n = 0 and n = — 1: This sequence of 


poles is regularly distributed on the complex time 
plane. Furthermore it can be shown that a polar 
singularity corresponds to each local maximum (or 
minimum) of the real time solution x p (t). Such 
regular distribution of singularities reflects the 
corresponding periodic behavior of the real time 
solution. What happens to the distribution of the 
singularities when the parameters o, 6, c of the Duff¬ 
ing equation are not zero? If these parameters are 
chosen such that the behavior of its solution x(z) on 
the real time axis is chaotic, then the corresponding 
sequence of singularities in the complex time plane 
associated to the local extrema becomes very irreg¬ 
ular. A complete study on the properties of the 
solutions of the Duffing equation on the complex 
time plane can be found in [Konno &; Tateno, 1984]. 
Other numerical investigations were performed also 
on the Lorenz equations [Tabor & Weiss, 1981]. 
In the case of the periodic regime of the Lorenz 
equations (limit cycle) the arrangement of singu¬ 
larities (poles) reflects, as for the Duffing equa¬ 
tion, the corresponding periodicity of the real time 
solution. As the dynamical regime goes toward the 
chaotic one the corresponding arrangement of sin¬ 
gularities becomes very irregular. Moreover, as for 
the Duffing equation, it was numerically shown that 
a singularity (the nearest to the real time axis) can 
be associated to each position t on the real time 
axis where a local maximum (or minimum) of the 
real time solution occurs. It was shown that the 
distances of these singularities from the real time 
axis are related to the real values of the solution 
[Tabor &: Weiss, 1981; Konno &; Tateno, 1984]. 
In other words, let be r = Imz the imaginary 
part of the position z = t + e + it (0 < |e| <C 1) 
of the singularity corresponding to a local maxi¬ 
mum (or minimum) of /(f), then as |r| decreases 
the value of |/(f)| increases (and viceversa). From 
these last remarks it follows that the sequence of 
points {tj, each associated to a local max¬ 

imum (or minimum) of the function /(t), can be 
thought of as a representation of the sequence 
of singularities of /(z) nearest to the real time 
axis. A point to be stressed is that for a sys¬ 
tem of ODEs the position of the singularities (if 
they exist) is deterministically defined [Ramani 
et al., 1989]. In other words the pattern of sin¬ 
gularities in the complex time plane is directly 
determined by the properties of the system of 
ODEs. For example, as seen before, for the ODE 
x + x 3 = 0 the set of singularities is well defined: 
These are poles and are disposed in a regular 
lattice on the complex time plane. These theoretical 
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and numerical results suggest the first algorithm, 
whose protocol is the following: 

(a) Given the time series s(t ) let us consider the 
set of couples {(tj,s(tj) : j = 1,n)} where s(tj) 
are the local maxima (or indifferently the local 
minima); 

(b) define the quantity 


L = 


n —1 

d 
1 


Sfe+i “ + [ s (^+i) - s fe)] 2 (!) 


Nj- 


that represents the length of the broken line 
joining the points {tj, s(tj ) : j = l,n}; 

(c) test if the value of L can be statistically 
accounted for by a linear Gaussian stochastic 
process obtained from the original time series 
s(t) by using Fourier Shuffled surrogate data. 
It is worth recalling that Fourier Shuffled surro¬ 
gates preserve the histogram of amplitudes and, 
approximately, the linear time correlations of 
the original data. More details on the surrogate 
data properties can be found in [Theiler et al., 
1992]. In what follows the above algorithm will 
be denoted as PSC (Pattern of Singularities in 
the Complex time plane). 


2.2. NET algorithm 


[Abarbanel et al., 1993]. Given a scalar time series 
s(t) - M(x), where x = F(x) and M represents 
the measurement function, the reconstruction the¬ 
orem [Takens, 1981] ensures that for a sufficiently 
high embedding dimension E the dynamical prop¬ 
erties of the original system are equivalently de¬ 
scribed by y = G(y), where y = [s(t), ds(t)/dt, ..., 
d E ~ 1 s(t)/d E ~ 1 t], y G R e , E > 2m + 1 (m is the 
box-counting dimension of the attractor). We can 
therefore concentrate our attention on the dynam¬ 
ical system y = G(y), y € R E . In principle the 
probability density function p{y,t) can be deter¬ 
mined and, from it [Soong, 1973], the joint proba¬ 
bility density p(s,s,'s,t) can be obtained (here we 
are supposing E > 3). Now, according to [Soong, 
1973] the expected number of extrema for unit time 
of the time series is given by 

/ +oo r+ oo 

ds / p(s,0,s,t)ds. (2) 

-oo J —OO 


In the particular case of a stationary process, 
q(t) is independent of t. For example, if the proba¬ 
bility density is Gaussian with zero mean it can be 
easily shown that 


. . 1 cr(s) 

9<!) = *> = 


(3) 


Several arguments support the correctness of a 
statistical description of a chaotic dynamical 
system x = F(x), x € R n that, as it is well 
known, exhibits sensitivity to the initial conditions 
[Eckmann & Ruelle, 1985]. Among them the most 
important is that in nature the measurement pro¬ 
cess, by which the observer interacts with a physical 
system, is characterized by a finite precision. Con¬ 
sequently, the state of the system at time t is not 
represented by a point in phase space but, rather, by 
a small region whose size reflects the finite precision 
of the measurement. Other sources of delocalization 
of the state are the incomplete specification of the 
initial conditions or the roundoff errors in numerical 
calculations. For a chaotic dynamical system, rep¬ 
resented by x = F(x), the probability density 
function p(x,t ) of a given system state x can be 
determined by solving the Liouville equation which 
is a particular case of the differential Chapman- 
Kolmogorov equation describing Markov processes 
[Gardiner, 1983; Nicolis, 1995; Soong, 1973]. If 
we have at our disposal only a time series, the 
probability density can be evaluated numerically 


where o(s) and a(s) are, respectively, the stan¬ 
dard deviation of s and s. By assuming that the 
dynamical system is ergodic we can estimate 
q(t) from the experimental time series without 
knowing p(s,s,s,t). We remark that in Eq. (2) 
p(s, 0, s, t)dsds represents the probability to find the 
intersection of the trajectory, described by y, within 
the rectangle of sides s + ds and s + ds in the plane 
(s, s) of the Poincare section defined by s = 0. 
If the evolution law of y is characterized only by 
linear time correlations we expect that the statistic 
q{t), computed using surrogata data preserving 
linear time correlations, will be not statistically 
different from that obtained with the original 
data. 

Consequently the idea is to test whether the 
number of extrema for unit time of an experi¬ 
mental time series s(t) can be accounted for by a 
linear Gaussian stochastic process. If this does not 
happen then we have an indication that the gen¬ 
erating process responsible for the observed s(t) is 
nonlinear [Theiler et al., 1992], This suggests that 
the probability density function p(s,s,s,t ) cannot 
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be associated to a linear stochastic process. The 

recipe of the algorithm is the following: 

(a) The number of extrema for unit time vq of the 
given time series s(tj) (j = 1, ..., N) is deter¬ 
mined and used as discriminating statistic; 

(b) from the time series s(tj) (j = 1 , ... , N) a set 
of n surrogate time series s s j (t) is generated; 

(c) the average number of extrema for unit time 
of the surrogate time series v s = ^ YHj =i u s,j is 
computed 

(d) by using the standard deviation a u of the set 
u St j (j = 1,2, ..., n), i^o and u 3 are compared 
to test if they are statistically different [Theiler 
et al., 1992]. Hereupon to indicate this al¬ 
gorithm we will use the abbreviation NET 
(Number of Extrema for unit Time). 


3. Numerical Results 

A number of time series, from well-known models, 
were generated to test the PSC and NET algorithms 
and to find their possible limitations. In the pre¬ 
vious section we considered only continuous time 
dynamical system. A naturally arising question is 
whether the NET and PSC algorithms work also 
for maps. There is not a clear answer to this ques¬ 
tion from the theoretical point of view. Therefore 
we decided to investigate numerically the perfor¬ 
mance of both algorithms also on the well-known 
Henon map. 




♦ lorenz 
■ henon 

▲ lorenz+nolse 
x henon+nolse 
X slnewave 

• gaussian noise 
+ coloured noise 
- slnewave+noise 


The Lorenz equations, 
x = cr(y — x) 

i, = px - y - xy + e{(t) (4) 

i = -j3z + xy 

with parameters a = 10, p = 28, (5 = 8/3 corre¬ 
sponding to the chaotic regime, were integrated nu¬ 
merically; in Eq. (4) £( t ) is the Gaussian white noise 
of unit standard deviation. Two time series were 
obtained for the case e = 0 and e / 0, respectively. 
Two other time series were generated by using the 
Henon map: one without noise and the other with 
additive noise. Two time series, representing Gaus¬ 
sian white noise and colored noise, were also gener¬ 
ated. Finally the last two time series were obtained 
from the signal s(t) = sin (uit) + e£(f) in the cases 
e = 0 and e / 0. The occurrence of a local maxi¬ 
mum (or minimum) in the time series s(t ) was deter¬ 
mined from the change of sign of the first difference 
series. In Fig. 1 the results obtained using the NET 
algorithm are plotted in terms of the significance 
S = \vq — v s \/&v [Theiler et al, 1992]. For the left 
panel N = 1000 points were used and N = 2000 
for the right one. For both panels 10 phase ran¬ 
domized surrogate data were used for each time se¬ 
ries. According to [Theiler et al., 1992] for S values 
> 2 the null hypothesis that the data came from a 
linear Gaussian stochastic process can be rejected 
with a probability greater or equal to 95%. In Fig. 1 



♦ lorenz 
■ henon 

A lorenz+noise 
X henon+nolse 
X sinewave 

• gaussian noise 
+ coloured noise 
- sinewave+noise 



Fig. 1. NET algorithm: Significance S = \i>o — v s \ /ov for the model time series. In what follows £(t) is Gaussian white noise 
of unit standard deviation. The noisy Lorenz time series was obtained with e = 0.2; for the sinus wave sin(wt) is lj = 0.1; for 
the sinus wave plus noise sin(t ot) + e£(f) is w = 0.1 and e = 0.3; the colored noise was obtained from x(t + 1) = ax(t ) + e£(t) 
with a = 0.85 and e = 0.2; the amplitude of the additive noise in the Henon map was e = 0.36. The sampling time for the 
Lorenz equation is At = 0.025. (n = 10), (a) N = 1000, (b) N = 2000. 
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Fig. 2. As in Fig. 1 but using n — 20 surrogates for each time series, (a) N = 1000, (b) N = 2000. 
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Fig. 3. As in Fig. 1, but using Fourier Shuffled surrogate time series (n = 10), (a) N = 1000 (b) N = 2000. 


each symbol identifies a model time series as Qspeci- 
fied in the legend (this correspondence will be main¬ 
tained). The test gives a clear indication of nonlin¬ 
earity only for the time series Lorenz, Lorenz plus 
noise, Henon map and Henon map plus noise that 
are known to be nonlinear dynamical systems. The 
case of the sinus wave is more problematic: The S 
value is slightly above 2 for both panels and this 
disagrees with the fact that the time series is lin¬ 
ear. However, as will be shown in the following, 
this anomalous behavior is corrected when the PSC 
algorithm is applied. In Fig. 2 the results obtained 
with the NET algorithm, and using 20 surrogate 
data for each time series, are plotted (the num¬ 
ber of points used for each time series are as in 
Fig. 1). The results are qualitatively similar to 
those in Fig. 1. However, from the left panel of 
Fig. 2 the noisy Lorenz time series seems to be 


more nonlinear than the corresponding noise-free 
one. But, as the number of points is-increased (from 
1000 to 2000) this discrepancy disappears (right 
panel). The analysis was repeated using 10 sets 
of Fourier Shuffled surrogate data for each model 
time series and the corresponding results are plot¬ 
ted in Fig. 3. These last results are very similar 
to those obtained with the phase randomized sur¬ 
rogate data; this means that the computed statistic 
is insensitive to slight changes in the distribution of 
the amplitudes of the time series. Next the PSC 
algorithm was applied to the same model time se¬ 
ries. The PSC algorithm requires the calculation 
of the length L, defined in Eq. (1), on the origi¬ 
nal time series and the statistical evaluation of its 
difference from L s = ^ Y%= i Lj ( Lj is calculated 
on the jth surrogate time series). To identify the 
local extrema (maxima in this case) we used the 
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Fig. 4. PSC algorithm using local maxima (n = 10), (a) N = 1000, (b) N = 2000, S — \L — L s \/ctl- Fourier Shuffled 
surrogata data are used. 


same protocol as for the NET algorithm. In Fig. 4 
the values of the significance S = \L — L s \/ai ob¬ 
tained using the PSC algorithm are plotted (ol is 
the standard deviation of the set {Lj}). The num¬ 
ber of points in the left panel is N = 1000 and 
N = 2000 in the right one. The time series corre¬ 
sponding to the sinus wave plus noise seems to show 
a nonlinear character (left panel). However, as the 
number of points increases this apparent nonlinear¬ 
ity disappears (right panel). Figure 5 shows the 
results obtained with the PSC algorithm by using 
the local minima instead of maxima. 

The results obtained in the Henon map, Henon 
map plus noise and Gaussian white noise increas¬ 
ing the series length up to IV = 4000 points are 
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Fig. 5. PSC algorithm: The local minima are used (n = 10), 
N = 2000, S = \L — L s \/<jl. Fourier Shuffled surrogata data 
are used. 


reported in Table 1. Similar results were also ob¬ 
tained for the other time series (data not shown). 
The conclusion arising from these numerical exper¬ 
iments is that increasing the length of the time 
series does not produce results contrasting with 
those obtained using shorter sequences of data. Fur¬ 
thermore in presence of nonlinearity (S values > 2) 
the value of S improves (statistically) as the length 
of the time series increases. 

Another important point that we faced is 
whether the PSC and NET algorithms were able to 
detect nonlinear correlations in time series coming 
from a high dimensional chaotic attractor. Some 
preliminary considerations are needed before show¬ 
ing the numerical results obtained with both meth¬ 
ods. As shown in the previous section the PSC 
algorithm requires the computation of L defined by 
Eq. (1). The set of points : j = l,n)}, 

where s(tj ) are the local maxima (or minima) of the 
time series s(t), is a bidimensional Poincare section 
[Hegger & Kantz, 1997; Kantz & Schreiber, 1997] 
defined by s(t) = 0, s(t ) < 0 (or s(t) = 0, s(t) > 0 
for the minima). For high dimensional attractors 
(dimension > 3) the Poincare section can mimic the 
one due to a stochastic process [Kantz & Schreiber, 
1997]. Similar considerations are also valid for 
the NET algorithm. Therefore it is expected that 
both PSC and NET algorithms can have problems 
with time series from a high dimensional attractor. 
We performed numerical experiments by using time 
series generated with the Mackey-Glass 

i(f) = - ta(() (5) 
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Table 1. Values of significance S obtained with the PSC and NET algorithms against 
the number of points N of the time series. 


Number of points 

Henon Map 

Henon Map + Noise 

Gaussian White Noise 

N 

PSC 

NET 

PSC 

NET 

PSC 

NET 


3.6 

4.8 

2.7 

2.86 

0.7 

1.01 

2000 

6.4 

3.6 

2.1 

2.87 

0.1 

1.1 

3000 

6.6 

6 

2.9 

3.95 

0.08 

0.46 

4000 

13.4 

9.5 

4.3 

4.7 

0.49 

0.09 


time-delay differential equation [Mackey & Glass, 
1977]. In this study we set a = 0.2,6 = 0.1 ,c = 10. 
Two time series of N = 4000 points correspond¬ 
ing to r = 30 and r = 100 were generated; the 
sampling time was At = 10. The corresponding 
correlation dimension of the attractor is ~ 3 (for 
r = 30) and ~ 7.1 (for r = 100). The values 
of the significance S computed for these two time 
series are the following: (a) (r = 30): PSC algo¬ 
rithm S = 0.04, NET algorithm 5 = 2; (b) (r = 
100): PSC algorithm S = 0.02, ET algorithm S = 
0.6. The above results show that both algorithms 
have problems to detect nonlinear correlations in 
time series from high dimensional attractors. The 
NET algorithm, however, seems to be less sensi¬ 
tive to the high dimensional dynamics than the 
PSC one. 


4. Concluding Remarks 

The time series corresponding to chaotic dynamics 
(time series: Lorenz, Lorenz plus noise, Henon map 
and Henon map plus noise) exhibit the expected 
nonlinear character for both algorithms (5 > 2). 
Similarly the time series corresponding to Gaus¬ 
sian white noise and colored noise present values 
of the significance S below 2 for both algorithms. 
This means that the null hypothesis that they are 
generated by a linear stochastic process cannot be 
rejected with a probability higher than 95%. The 
time series corresponding to the sinus wave, as 
seen in Figs. 1-3, exhibits an unexpected nonlinear 
behavior using the NET algorithm. This behavior 
disappears for the sinus wave and appears for the 
sinus wave plus noise (left panel of Fig. 4) when the 
PSC algorithm is used. However, as the number of 
points is increased (right panel of Fig. 4) this effect 
for the sinus wave plus noise disappears. We studied 
the dependence of the values of the significance S 
from the length of the time series. The results were 


that the growth of the number of points N of the 
sequence produces values of S coherent with those 
obtained using shorter data sets; furthermore in the 
presence of nonlinear correlations in the time series 
(S values > 2) the value of S improves statistically 
as N increases. 

We checked also the performance of both algo¬ 
rithms to detect nonlinear correlations in time series 
from a high dimensional attractor. The result was 
that both algorithms fail in such a case; furthermore 
it seems that the NET algorithm is less sensitive 
than the PSC one to the growth of the dimension 
of the attractor. We remark that, using the classi¬ 
cal algorithm to estimate the correlation dimension 
from time series, the high dimensional dynamics of 
the Mackey-Glass system is detected, using large 
amount of data points (»4000), correctly [Ding 
et al. , 1993]. It is presently under investigation to 
assess whether some change in the protocols of both 
algorithms can improve their perdormance on time 
series from a high dimensional attractor. 

All the above results lead to the following 
remarks: (a) For low dimensional chaotic dynamics 
the nonlinear property (S > 2) is detected by both 
algorithms; (b) for colored noise or Gaussian white 
noise significance values S < 2 are obtained with 
both algorithms; (c) when the signal power spec¬ 
trum exhibits a strong sharp peak the two algo¬ 
rithms give contrasting results. This can be par¬ 
tially explained by the fact that the use of surrogate 
data for periodic signals with long coherence time 
can lead to ambiguous results that need to be inter¬ 
preted very carefully [Theiler et al., 1992]. These 
preliminary results suggest that the two methods 
used together are able to discriminate between 
low dimensional chaotic dynamics and noise. The 
last remark is that the NET and PSC algorithms 
do not require the reconstruction protocol of the 
phase space as it happens for many of the current 
tools for the nonlinear analysis of time series (Lya¬ 
punov exponents, dimension calculations, nonlinear 
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prediction methods, etc.) [Abarbanel et al. , 1993]; 
consequently we do not need to determine the 
embedding dimension E and the lag-time r. More¬ 
over these more complex tools to study time 
series are not trivial to be implemented using nu¬ 
merical codes. On the contrary the NET and PSC 
algorithms are very simple to implement and they 
do not require large amount of data in order to 
check the presence of nonlinear correlations in time 
series from low dimensional attractor. Therefore 
the NET and PSC algorithms can be used as a 
first test before using the aforementioned more 
complex tools. 
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Secure communication via chaotic synchronization is demonstrated using dynamical systems 
governed by delay deferential equations. Strange attractors of such systems can have an arbi¬ 
trarily large number of positive Lyapunov exponents giving rise to very complex time signals. 
This feature can provide high security of masked messages. 


1. Introduction 

Synchronization of chaotic systems has aroused 
much interest in recent years, particularly in light 
of potential application of this phenomenon in se¬ 
cure communication. Different approaches for con¬ 
structing synchronized chaotic systems are pro¬ 
posed. Among them is the approach of Pecora and 
Carroll [1990] who show that when a state variable 
from a chaotically evolving system is transmitted 
as an input to a replica of part of the original sys¬ 
tem, the replica subsystem (receiver) can synchro¬ 
nize to the original system (sender). To use this 
phenomenon for masking messages one can trans¬ 
mit to the receiver a summary signal consisting of 
a chaotic signal generated by sender and a small 
information signal containing a message [Kocarev 
et al. , 1992], Then the message can be recovered 
by synchronizing the receiver with the scalar signal 
which is transmitted from the sender. The short¬ 
coming of this approach is that the information 
signal violates the synchrony between sender and 
receiver giving rise to reconstruction errors in the 
recovered information. 

Recently, Kocarev and Parlitz [1995] have pro¬ 
posed a generalization of the Pecora and Carroll 
[1990] method, which extends the capabilities for 
constructing synchronized chaotic systems. This 


new approach enables the information signal to 
be integrated in the sender as a driving signal. 
The scalar signal which is transmitted from the 
sender to the receiver is a function of the sender 
state variables and the information signal. At ideal 
conditions the information signal can be recovered 
exactly, without the reconstruction errors. 

Most theoretical as well as experimental stud¬ 
ies of secure communication reported so far concern 
low-dimensional systems with one positive Lya¬ 
punov exponent. Perez and Cerdeira [1995] have 
shown that messages masked by such simple chaotic 
processes, once intercepted, can be sometimes read¬ 
ily extracted. To improve security it is desirable 
to use high-dimensional systems with multiple pos¬ 
itive Lyapunov exponents (hyperchaotic systems) 
which take advantage of the increased randomness 
and unpredictability. Recently, Peng et al. [1996] 
have demonstrated the possibility of synchronizing 
hyperchaotic systems by transmitting just a single 
scalar variable composed of a linear combination of 
state variables of the sender. Another way of syn¬ 
thesizing synchronized hyperchaotic systems using 
a series of low-dimensional subsystems as building 
blocks has been proposed by Kocarev and Parlitz 
[1995]. These methods have many advantages, but 
they have been applied only to finite-dimensional 
systems described by ordinary differential equations 
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(ODE’s). The number of positive Lyapunov expo¬ 
nents of such systems is limited by dimension of the 
state space. 

An alternative efficient approach of construct¬ 
ing synchronized hyperchaotic systems can be based 
on delay differential equations (DDE’s). Systems 
governed by DDE’s have an infinite-dimensional 
state space and can produce hyperchaos with an 
arbitrarily large number of positive Lyapunov ex¬ 
ponents. A typical example of these is the Mackey 
and Glass (MG) [1977] system: 

x = f[x(t - r),a] - cx, (1) 

where f[x(t — r),a ] = ax(t — r)/[ 1 + x b (t — r)] 
and a = {a, 6} denotes a set of parameters of the 
function /. All parameters a, b , c, and r are sup¬ 
posed to be positive. We choose this model for our 
studies since it has been thoroughly investigated in 
the literature (e.g. [Farmer, 1982]) and is easily im- 
plementable electronically [Namajunas et al., 1995]. 
The number of positive Lyapunov exponents as well 
as dimension of strange attractor of this system can 
easily be controlled by varying the delay time r. 
For fixed values of the parameters a, b, and c, both 
these quantities increase linearly with the increase 
of r [Farmer, 1982]. 

Below we demonstrate secure communication 
with two types of senders based on a single MG 
system and using a set of coupled MG systems. 

2. Communication Based on a Single 
MG System 

Let us consider the communication system with the 
equations of the sender, the transmitted signal, and 
the receiver given by 

x — f[x(t — t) + i(t), a] — cx sender, (2) 

s(t ) = x(t — t) + i (t) transmitted signal, (3) 

V = /[ s (£)>«] - cy receiver , (4) 

where i(t ) is the information signal. The sender (2) 
is an infinite-dimensional system described by non¬ 
linear DDE, while the receiver (4) is presented by 
a simple linear ODE driven (through a nonlinear 
function /) by transmitted signal s(t). Subtracting 
Eq. (4) from Eq. (2) we obtain for the difference 
e = x — y a simple linear ODE e = — ce. This equa¬ 
tion possesses an unique globally stable fixed point 
e = 0. Thus the synchronized state x — y is globally 
stable, i.e., the receiver synchronizes to the sender 


independently on initial conditions. The character¬ 
istic time of synchronization is given by 1/c. At the 
receiver the information iR can be recovered as 

iR = s(t) - y{t - t ) . (5) 

Figure 1 illustrates secure communication for 
the parameter values of the MG system a = 0.2, 
b = 10, c = 0.1, and t = 100. Here and below 
we perform numerical integration of the underlying 
DDE’s by a second order Runge- Kutta method tak¬ 
ing an integration step h = 0.1. Without the infor¬ 
mation signal, i(t ) = 0, the sender [Eq. (2)] has five 
positive Lyapunov exponents. The eleven largest 
exponents multiplied by factor 10 3 are {3.53, 2.88, 
2.09, 1.28, 0.46, 0.00, -0.54, -1.49, -2.63, -3.83, 
—4.98). This corresponds to information dimen¬ 
sion [Kaplan & Yorke, 1979] of the strange attractor 
equal to d\ 10.35. 

Figure 1(a) shows the transmitted signal s(t) in 
the case of a sine information signal i(t) with the 
period T = 80 and amplitude A = 0.1 switched on 
at the moment t = 200 and switched off at t = 800. 
The recovered information signal iR{t) is presented 
in Fig. 1(b). It coincides with the original infor¬ 
mation signal i(t), to within the error of numerical 
integration. 

In a real experiment, the parameters of sender 
and receiver are not exactly the same. This factor 
and additional noise in a transmission channel will 
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t 

Fig. 1. Numerical simulation of a communication scheme 
based on a single MG system in the case of sine information 
signal i = 0.1sin[27r(t — 200)/80] switched on at the mo¬ 
ment t = 200 and switched off at t = 800. The parameter 
values of the MG system are a = 0.2, b = 10, c = 0.1, and 
r = 100. (a) Transmitted signal s. (b) Recovered information 
signal iR. 
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result in reconstruction errors. However, due to the 
global stability of the synchronized state x = y, 
the system is robust with respect to these factors. 
Suppose that the parameters of the driving and the 
response systems are different, say a and a', re¬ 
spectively. If the difference 5a — a — a' is small, 
the difference e = x — y is governed by a linear 
nonhomogenous differential equation 


e — 5a 


dj_ 

da 


— ce. 


Clearly the difference e is proportional to the pa¬ 
rameter mismatch 5a, e oc 5a. The same is valid 
for the difference between the original and the re¬ 
constructed information signal, i — in oc 5a, which 
defines the reconstruction errors. The information 
signal can be properly recovered if the difference 
5a is small compared to the amplitude of informa¬ 
tion signal. The influence of the mismatch between 
parameters to the reconstruction errors is numeri¬ 
cally illustrated in Fig. 2(a). The parameter a is 
taken to have different values at the sender and the 
receiver, a and a 1 , respectively. We measure the 
re constructio n errors by the r.m.s. deviation o = 
y/({i — ir) 2 ), where (•) denotes the time average. 
As is seen from the figure, the reconstruction er¬ 
rors decrease linearly with the decrease of the dif¬ 
ference 5a = a — a'. Similar results are observed 
in the presence of noise. To model an influence of 
noise in the transmission channel, we add (at ev¬ 
ery step of numerical integration) random numbers 
uniformly distributed in the interval [— a n /2,a n /2] 
to the transmitted signal s(t). Figure 2(b) shows 
that for small noise the reconstruction errors are 
proportional to the noise amplitude a n . 
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Fig. 2. Dependence of the r.m.s. deviation a = 
V<(i~ *«) 2 > on (a) parameter mismatch 5a and (b) noise 
amplitude a n ■ The parameters of the MG system are the 
same as in Fig. 1. The information signal represents the sine 
wave i = 0.1 sin(27rt/80). 


Although the MG system allows us to gener¬ 
ate very complex signals, it is still described by 
a very simple DDE. In the absence of informa¬ 
tion signal, i(t ) = 0, such simple DDE can be 
reconstructed from transmitted signal s(t ) using, 
for example, the method proposed by Biinner 
et al. [1996]. Certainly, the presence of information 
signal (i ^ 0) considerably complicates the possibil¬ 
ity of such reconstruction, since the sender becomes 
a nonautonomous system. To make the communi¬ 
cation system even more secure one can increase the 
complexity of the sender by using a set of coupled 
MG systems. 


3. Communication Based on a Set of 
Coupled MG Systems 


To illustrate the possibility of more complex con¬ 
structions of a communication system we consider 
the specific example with the equations of the 
sender, the transmitted signal, and the receiver 
given by 


{ &l = f[xi(t~n) + i(t),a 1 ]-cixi 

Sa.ux{t) = Xl(t-Ti)+ i(t) 

X 2 =f[x 2 (t-T 2 ) + ks anx (t),a 2 ]~C 2 X 2 


> sender, 


( 6 ) 


s(t ) = x 2 {t — t 2 ) + ks aux (t ) transmitted signal. 


S /2 = f[s(t),a 2 ] - c 2 y 2 
Saux(i) = [s{t) - y 2 (t - r 2 )\/k 
in = f[sa.ux(t),ai] - CiyI 


(7) 

> receiver. (8) 


The idea [Kocarev &: Parlitz, 1995] behind this con¬ 
struction is based on successive transmission of in¬ 
formation signal through different building blocks 
(here they represent the MG systems) of the sender 
and then successive recover of the transmitted sig¬ 
nal at the different blocks of the receiver. Here we 
consider the construction with only two building 
blocks. The generalization for an arbitrary num¬ 
ber of building blocks is straightforward. 

This construction as well as previous examples 
provides the global stability of the synchronized 
state x\ = yi, x 2 = y 2 . Indeed, the difference 
e 2 = x 2 — ?/2 is governed by the linear equation 
e 2 = —C 2 e 2 possessing the globally stable fixed 
point e 2 = 0. This guarantees the global stabil¬ 
ity of the state x 2 = y 2 . Taking into account the 
last equality we obtain s aux (t) = s aux (t), and the 
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equation for the difference ei = X\ — y\ takes the 
form ei = —C\e\. Thus the state x\ = yi is also 
globally stable. At the receiver the information i R 
can be recovered as 

i R = [s(t) - y 2 {t - r 2 )\/k - yi(t - n). (9) 

Figure 3(a) shows the amplitude spectrum of 
the transmitted signal in the case of sine informa¬ 
tion signal i(t ) = 0.1 sin(27ri/40). The amplitude 
spectrum of the recovered information signal i R is 
shown in Fig. 3(b). The information signal is prac¬ 
tically invisible in the spectrum of the transmitted 
signal. 

4. Conclusions 

We have demonstrated an efficient way of construct¬ 
ing secure communication systems using delay dif¬ 
ferential equations. The sender of such systems has 
an infinite-dimensional state space and can produce 
chaotic carrier signals with an easily controllable 



Frequency (Hz) 


Fig. 3. Numerical simulation of a communication scheme 
based on a set of coupled MG systems in the case of sine 
information signal i = 0.1 sin(27rf/40). The parameter values 
of the MG systems are ai = 0.22, b\ = 12, ci = 0.1, ri = 110, 
d 2 = 0.2, 62 = 10, C 2 = 0 . 1 , 72 = 100, and k = 0 . 1 . (a) Am¬ 
plitude spectrum of the transmitted signal s. (b) Amplitude 
spectrum of the recovered sine signal iR. 


number of positive Lyapunov exponents. The re¬ 
ceiver is constructed as a linear low-dimensional 
system. The synchronization with the sender is 
achieved by transmitting only single scalar signal. 
An important advantage of this construction is the 
global stability of the synchronized state. As a 
result these communication systems are insensi¬ 
tive to small noise levels and mismatching between 
parameters. 
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An experiment on Benard-Marangoni time-dependent convection shows evidence of an am¬ 
plitude turbulent regime in the temperature signal which is modeled by a delayed dynamical 
system. Application of a control procedure, which perturbs the value of the delay time, leads 
to the control of such dynamical regime, by suppression of phase defects and stabilization of 
the regular oscillations. The control technique is robust against the presence of large amounts 
of noise. 


1. Introduction 

The original idea of Ott, Grebogi and Yorke [Ott 
et al., 1990] on chaos control has generated many 
different theoretical schemes and experimental ap¬ 
plications facing the problem of controlling unstable 
periodic orbits (UPO’s) in chaotic concentrated sys¬ 
tems (CS), that is in systems modeled by ordinary 
differential equations. 

Controlling spatially extended systems (ES), 
i.e. systems ruled by partial differential equations 
whose order parameter y is a m-dimensional vec¬ 
tor (m > 1) in phase space, with k components 
(k > 1) in real space, is still an open problem. Even 
though some proposals have been put forward for 
the case k = 2 [Lu et al. , 1996], experimentally im- 
plementable tools have not yet been introduced for 
the control of unstable periodic patterns (UPP) in 
extended systems. 

The link between CS and ES is provided 
by delayed dynamical systems (DS), i.e. systems 


ruled by 

y = ?{y, ya ), (i) 

y = y(t), dot denotes temporal derivative, F is a 
nonlinear function, yd = y(t — T), and T is a time 
delay. 

The evidence of the analogy between DS and ES 
was given experimentally for a CO 2 laser with de¬ 
layed feedback [Arecchi et al., 1992] and supported 
by a theoretical model [Giacomelli & Politi, 1996]. 

The DS to ES conversion is based on a two vari¬ 
able time representation, defined by t — a + 6T , 
where 0 < cr < T is a continuous space-like variable 
and 6 is a discrete temporal variable [Arecchi et al., 
1992], In this framework, the long range interac¬ 
tions introduced by the delay can be reinterpreted 
as short range interactions along the 0 direction 
(yd = y{cr, 6 — 1)) and the formation and propaga¬ 
tion of space-time structures, as defects and/or spa- 
tiotemporal intermittency can be identified [Arecchi 
et al., 1992; Giacomelli & Politi, 1996]. 
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For delays T larger than the period of oscilla¬ 
tion of the system, the behavior of a DS is analo¬ 
gous to that of an ES with k — 1. Namely, DS may 
display phase defects, i.e. points where the phase 
suddenly changes its value and the amplitude goes 
to zero. In this paper we show evidence of these 
phase defects in a recent experiment on Benard- 
Marangoni convection [Mancini & Maza, 1997], and 
we propose a control technique to suppress them. 
The control restores regular patterns within an 
amplitude turbulent regime, which implies the pres¬ 
ence of a large number of defects. The control effi¬ 
ciency persists even in presence of a large amount 
of noise. 


2. The Experiment and 

The Delayed Dynamical Model 

The experimental setup is depicted in Fig. 1. A 
cylindrical cell (diameter 128 mm) confines a fluid 
layer of silicon oil (Pr ~ 3000) with the free sur¬ 
face open to the atmosphere and heated from the 
bottom. The heater does not cover the whole of the 
container giving open boundaries for the heating. A 
convective instability driven by buoyancy and tem¬ 
perature dependent surface tension (80% and 20%, 
respectively), takes place as the heating is increased 
giving rise to a stationary planform [Fig. 1(a)]. This 
pattern is composed of four convective cells located 
on the heater region, but the flow is developed over 
all the container size. 

Following an imaginary drop of fluid traveling 
with the flow in one of these cells, the drop is heated 
while traveling near the bottom over the heater, 
rises up to the centre, travels out near the surface 
until it becomes cold and then falls down near the 
lateral boundaries. Finally, the drop is fed back to 
the heater region completing a round trip in a mean 
time T [Fig. 1(b)]. 

If the heating is further increased a time- 
dependent regime arises, generating spatio- 
temporal modulations of the stationary velocity 
and temperature fields. The origin of this behavior 
is related with a thermal boundary layer instability 
which give rise to thermals or “hot plumes” which 
are dragged by the flow along the cell. This behav¬ 
ior can be seen in the space-time image of Fig. 1(c). 
An experimental measurement of the temperature 
at the center of the cell shows modulated oscilla¬ 
tions which have a power spectrum composed by 
two frequencies clearly differentiated (plus their 



Fig. 1. (a) Image of the stationary planform below the time- 

dependent regime, (b) Cross-section of the experimental 
setup. Above the time-dependent threshold, the thermals 
coming from the thermal boundary layer generate hot drops 
which are dragged by the flow along the convective cells in a 
mean time T. (c) Time-dependent regime observed from the 
planform. The white traces in the spatiotemporal diagram 
correspond to the effect of the hot drop traveling near the 
surface in the x-axis direction. See [Mancini & Maza, 1997] 
for futher details. 

nonlinear combinations terms). One of these 
frequencies corresponds to a relaxation oscillation 
inside the thermal boundary layer, the other one 
corresponds to the characteristic time necessary for 
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Fig. 2. (a) Experimental time behavior of the temperature signal. Vertical axis reports the temperature at the center of the 

cell near the bottom and horizontal axis reports the time in seconds (T ~ 330 sec.), (b) Expanded view of the signal within 
the box which exhibits a phase jump (indicated by the arrow). 


a round trip of a thermal along the convective 
cell. Further details on the experiment and a de¬ 
tailed discusion of this mechanism can be found in 
[Mancini k. Maza, 1997]. 

If the temperature of the heater is further in¬ 
creased, a chaotic regime is reached. In this regime, 
an experimental measurement of the temperature at 
the center of the cell yields the data in Fig. 2. The 
signal shows trains of modulated oscillations, inter¬ 
rupted by localized events (phase defects), wherein 
the phase changes suddenly and the amplitude de¬ 
creases to zero. Figure 2(b) highlights the presence 
of a phase defect within the data. 

The experimental configuration provides a 
natural delayed interaction between thermals and 
thermal boundary layer since it reiterates at each 
position the local value of the order parameter 
after a time delay T, which equals the time lag 
necessary for the trip of the cell. Moreover, it 
involves a pertubated state far from the time 
dependent convection threshold. We propose a 
nonlinear model for the experimental temperature 


signal: 

A = eA +Pi (jT° A 2 {t - A 

+ A 4 (t-t')f(t’)dt^A, (2) 

i = n (s - - kA 2 ^j . (3) 

Here, all quantities are real. A represents the 
temperature, £ is a time dependent control param¬ 
eter, Pi, P 2 , Pi, k are suitable fixed parameters, fi 
is a measure of the ratio between the characteristic 
time scales for A and e, and S is a measure of the 
power provided to the system. 

The relaxation oscillations of the temperature 
in Fig. 1(b) are represented by the normal form 
of a Hopf bifurcation [Eq. (2)], in which the satu¬ 
rating terms include a delayed function modulated 
by f(t f ) to account the delayed action of the ther¬ 
mal inside the convective cell. Equation (3) mod¬ 
els the slow evolution (fi < 1) of the linear gain e, 
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which is enhanced by the external heating (S) and 
depressed by the convective motion (— kA 2 ) which 
tends to uniformize top and bottom temperatures. 
In general f(t') is a Gaussian-type function which 
expresses the lateral heat diffusion of the thermal to 
the main flow. However, we will consider the case of 
a perfectly localized pulse using a Dirac delta func¬ 
tion located a t' = T. The model can be written 
as: 

A = £A + p 1 A 2 {t-T)A + /3 2 A 4 (t-T)A, (4) 

i = g ^S — — e — kA 2 ^j . (5) 

Even though Eqs. (4) and (5) have been here 
introduced for modeling a specific experimental 
situation (a chaotic transition associated with 
quasiperiodicity), they are in fact rather general. 
When T = 0, 5 < 0, pi > 0, & < 0, g > 0, 
pi > 0, k > 0, they model an excitable system, 
producing the so-called Leontovitch bifurcation, ev¬ 
idence of which has been shown experimentally on 
a CO 2 laser with intracavity saturable absorber 
[Plaza et al., 1997]. For T / 0, they are similar 
to the models already introduced to describe self- 
sustained oscillations of confined jets [Villermaux & 
Hopfinger, 1994], or memory induced low frequency 
oscillations in closed convection boxes [Villermaux, 
1995], or even the pulsed dynamics of a fountain 
[Villermaux, 1994], 

Adjusting pump and delay parameters (5 and 
T) in Eqs. (4) and (5), the system enters the chaotic 
region. This region, in fact, is split into two differ¬ 
ent regimes. For low T values, chaos is due to a 
local chaotic evolution of the phase, whereas no ap¬ 
preciable amplitude fluctuations are observed. This 
regime is called phase turbulence (PT). By increas¬ 
ing T, a transition toward amplitude turbulence 
(AT) is observed. In AT, the dynamics is dominated 
by the amplitude fluctuations, and a large number 
of defects is present. Both PT and AT have coun¬ 
terparts in a one-dimensional complex Ginzburg- 
Landau equation, for which the parameter space 
shows a transition from a regime of stable plane 
waves toward PT (Benjamin-Fair instability), fol¬ 
lowed by another transition to AT with evidence of 
space-time defects [Montagne et al, 1996] . 

3. The Control 

The aim of the present paper is to control AT by 
an adaptive technique recently introduced for chaos 
recognition [Arecchi et al., 1994], and applied to 


chaos control on CS [Boccaletti & Arecchi, 1995], 
chaos synchronization [Boccaletti et al, 1997], tar¬ 
geting of chaos [Boccaletti et al, 1997] and filtering 
of noise from chaotic data sets [Boccaletti et al, 
1997]. A direct application of such a technique 
to Eqs. (4) and (5) has been already provided in 
[Boccaletti et al, 1997]. In that case, a small con¬ 
tinuous perturbation U (t) of the local value of the 
temperature leads to the suppression of the phase 
defects, and restores the regular Hopf oscillations. 
Here we show an alternative strategy for the con¬ 
trol of AT, whereby tiny continuous modifications 
of the parameter T lead to a local control of the 
phase of the signal. In order to prove the efficacy of 
our method we use the system in Eqs. (4) and (5) 
in the AT regime. Here, the time delay is propor¬ 
tional to the spatial extension of the system. We 
will show that very small perturbations of the time 
delay are sufficient for the control of phase defects. 

Let us, therefore, consider the modified system 

A = eA + PiA 2 (t - (T + U(t)))A 

+ 02A 4 (t — (T + U(t)))A, (6) 

i = g — —e — kA 2 ^ . (7) 

The control algorithm which selects U(T ) can 
be summarized as follows. At time t n +1 = t n + r n 
(r n being an adaptive observation time interval to 
be later specified), the observer defines the variation 
A(t n+ 1 — Th ) — A{t n +\) between the actual and the 
delayed values of A (Th being the Hopf period). 
The corresponding variation rate 


An+i — — log 


A(t n+ 1 - Th) ~ A(t n+ 1 ) 
A(t n - Th) - A(t n ) 


( 8 ) 


allows to select a new time interval 


r n +1 = r n { 1 - tgh(fifA n+1 )), g > 0 (9) 

and consequently a new observation at the time 
t n+ 2 = t n+ 1 + r n+ i. In the following we perturb T 
by adding iteratively to it a controlling term given 
by 

U(t) = — (A(t-T H )-A(t)). (10) 

T n +1 

The details of the algorithm have been given in 
[Arecchi et al, 1994; Boccaletti et al, 1997], For 
relatively small perturbations, the following approx¬ 
imation holds. Let (r) denote the average of the 
{r n } set, then Eq. (9) can be written as 

Tn+ 1 - (r)(l - gXn+l) (11) 
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Fig. 3. (a) Space(<7)-time(0) representation of the controlling process for Eqs. (6) and (7). fi\ — 1, = —1/16, p = 0.8, 

Pi — 0.8, k = 11, S = 7, Th — 1.95. T = 50, AT regime. The dynamics is dominated by amplitude fluctuations, with the 
presence of defects. Phase defects appear as dislocations in such a representation. The algorithm (Ki = 0.3, K 2 = 0.07) 
suppresses the defects and restores the regular oscillation. Arrow indicates the instant at which control is switched on. 
(b) The behavior of the time delay T +U{t). 


a 



Fig. 4. T = 50, AT with 10% noise. Control with K\ = 0.3, 
K 2 = 0.07. Same stipulations and parameters as in the cap¬ 
tion of Fig. 2. Arrows indicate the instant at which control 
is switched on. 


where (i) r n has been replaced with its ensemble av¬ 
erage, and (ii) the tgh function has been linearized. 
In the same way, Eq. (8) can also be linearized as 


\(A - 

{ ) ~ (r) A(t) 


A(t - T h ) 
A(t - T H ) 


where the discretized stroboscopic observations 
have been approximated with a continuous inspec¬ 
tion. Combining Eqs. (11) and (12) into Eq. (10), 
this reduces to 


U(t) = K 1 (A(t-T H )-A(t)) 

+ K 2 (A(t-T H )-A(t)) (13) 

with K\ = l/(r) and K 2 = gj{r) 2 . The conse¬ 
quences of this approximation are relevant. First 
of all, for K 2 = 0 one recovers the Pyragas control 
method [Pyragas, 1992]. However, in our case, K\ 
and K 2 can be independently selected, and this in¬ 
troduces an extra degree of freedom with respect to 
[Pyragas, 1992]. 

In Fig. 3 we report the application of our 
method to Eqs. (6) and (7). The desired oscillation, 
which in the space-time representation gives rise to 
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a roll set, is controlled in AT [Fig. 2(a)]. Figure 2(b) 
reports the behavior of U(t), which is vanishing up 
to the instant at which control is switched on. 

Let us now discuss the robustness of our pro¬ 
cedure against external noise. For this purpose, we 
add white noise to the measured A values before 
the onset of the adaptive feedback control. In such 
a case, the noise does not act additively, since it af¬ 
fects the calculation of U(t), hence the local value 
of the time delay. As a consequence, the noise acts 
dynamically on the evolution of the system. A rel¬ 
evant result is that our method is robust against 
large amounts of noise, as it can be appreciated in 
Fig. 4 where the control is achieved within AT for 
10% noise. 

4. Conclusion 

We show that is possible to control delayed dynam¬ 
ical systems applying small perturbations on the 
time delay variable. The control algorithm is eas¬ 
ily implementable. The robustness of this method 
against noise has been verified. The proposed pro¬ 
cedure would imply an experimental setup wherein 
the control is achieved by modifying the cell length. 
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The stabilization of periodic solutions in the regime of spatiotemporal chaos through a diffusion 
parameter control is studied. We show that unstable plane waves in the Complex Ginzburg- 
Landau equation can be effectively stabilized in chaotic regimes such as phase turbulence and 
spatiotemporal intermittency or defect turbulence. 


1. Introduction 

There has been recently a considerable experimen¬ 
tal and theoretical effort to characterize Spatiotem¬ 
poral Chaos (STC) [Cross & Hohenberg, 1994; 
Gollub, 1994], Weak STC seems to be an ubiquitous 
phenomenon in large nonequilibrium systems. In 
some cases STC arise in the proximity of threshold 
and can be described within the context of weakly 
nonlinear theories. These theories are well devel¬ 
oped in the form of the so-called complex Ginzburg- 
Landau equations (CGLE) [Cross &; Hohenberg, 
1993]. The CGLE is a prototypical equation for 
a complex field A that exhibit STC [Chate, 1995]. 
It accounts for the slow modulations, in space and 
time of the oscillatory state in a physical system 
which undergoes a Hopf bifurcation [van Saarloos 
&; Hohenberg, 1992]. 


The control of spatiotemporal chaos is a com¬ 
plicated problem, and so, there is a wide variety 
of methods intended to control such chaotic behav¬ 
ior. There have been several attempts to achieve 
such control in the CGLE [Aranson et al., 1994; 
Bleich & Socolar, 1996; Mertens et al., 1994; 
Battogtokh & Mikhailov, 1996; Montagne & Colet, 
1997]. The most common approach is adding time- 
delayed feedback terms to the CGLE. The feedback 
can be either local [Bleich & Socolar, 1996] (at each 
spatial point, the field at the same point at pre¬ 
vious times is fed back) or global [Mertens et al., 
1994; Battogtokh & Mikhailov, 1996] (at each spa¬ 
tial point a term proportional to the integral of the 
field over the Spatial variable is fed back). Feedback 
has also been used for control in a nonlinear drift- 
wave equation driven by a sinusoidal wave [Gang, 
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1993]. And, in conjunction with a spatial filter, it 
has been applied to stabilize rolls and hexagonal 
structures in a model for a transversally extended 
three level laser [Lu et al., 1996] and to control 
filamentation in a model for wide aperture semicon¬ 
ductor lasers based on the Swift-Hohenberg equa¬ 
tion [Bleich et al., 1997]. 

A different approach using nonlinear diffusion 
effects has been introduced recently in [Montagne 
& Colet, 1997]. It was shown that an imaginary 
nonlinear diffusion term was able to stabilize unsta¬ 
ble plane waves in several regimes of STC. Here we 
generalize this study considering the case of having 
both, real and imaginary nonlinear diffusion terms 
in the CGLE. One of the main advantages of this 
generalization is the possibility of controlling plane 
waves outside the region where control works us¬ 
ing only an imaginary nonlinear diffusion term. In 
particular, in that case it was practically impossi¬ 
ble to stabilize plane waves with a large wavevec- 
tor, whereas this can be easily achieved in the case 
shown here. Also, in situations where control was 
already possible, this extension allows to achieve the 
stabilization of plane waves using a control parame¬ 
ter with smaller absolute value. As it happens with 
most of the control techniques, stabilization is pos¬ 
sible only when starting from an initial state close 
to the unstable orbit. If the initial condition is an 
arbitrary chaotic state, the system will explore the 
different regions in phase space of the chaotic at¬ 
tractor and eventually may approach the unstable 
orbit one wishes to stabilize, although this can take 
an extremely long time. If the system goes close 
enough to the unstable orbit then the control tech¬ 
nique shown here should work. Although we are 
not exploring in detail how close to the unstable or¬ 
bit the system has to be to achieve stabilization, we 
show numerically that our scheme is robust to finite 
size perturbations of the initial condition. 

In Sec. 2 we briefly describe the parameter 
regions for which different chaotic behaviors have 
been found for the CGLE and we introduce the 
modified equation. Section 3 is devoted to the lin¬ 
ear stability analysis of the plane wave solutions. 
We calculate for which parameter values the added 
term is able to stabilize plane waves in the STC re¬ 
gions of the CGLE. In Sec. 4 we show, integrating 
the equations numerically, that the region of stabil¬ 
ity of the plane waves when finite size perturbations 
are applied is in excellent agreement with the an¬ 
alytical prediction of the linear stability analysis. 
Finally we give some concluding remarks in Sec. 5. 


2. Model 

The one-dimensional CGLE [Cross & Hohenberg, 
1993; van Saarloos & Hohenberg, 1992; van Saar- 
loos, 1995; Newell et al., 1993] for a complex field 
A{x, t), describes the slow dynamics of spatially ex¬ 
tended systems close to a Hopf bifurcation, 

d t A = A + (1 + ia)d%A - (1 + ic 2 )\A\ 2 A . ( 1 ) 


We will assume periodic boundary conditions 
throughout the paper. This equation admits among 
other exact solutions, plane waves of the form 

A pw (x, t) = A 0 e^ kx -^ , (2) 

with amplitude Aq = \/l — k 2 , |A:| < 1 and fre¬ 
quency U) = C 2 + (ci — c 2 )k 2 . 

For 1 + C 1 C 2 > 0 plane wave solutions are lin¬ 
early stable for wave numbers smaller than a limit 
value |fc| < k,£ given by 

2 _ 1 + CiC 2 

E 3 + C\C 2 + 2 c\ 

For |fc| > kjg, plane waves are unstable to 
phase perturbations (Eckhaus instability [Eckhaus, 
1965]). The stability range vanishes at l + cic 2 = 0, 
the Benjamin-Feir-Newell (BFN) line, and there 
are no stable plane wave solutions for 1 + cic 2 < 0 . 

Numerical work for large system size [Shraiman 
et al., 1992; Chate, 1994, 1995] has identified re¬ 
gions of the parameter space displaying different 
kinds of regular and spatiotemporal chaotic behav¬ 
ior, leading to a “phase diagram” for the CGLE 
in the plane ci~c 2 . The five different regions, 
each leading to a different asymptotic phase, are 
shown in Fig. 1. Two of these regions are in the 
BFN stable zone and the other three in the BFN 
unstable one. In this paper we will concentrate 
in the regimes with a chaotic behavior, namely 
SpatioTemporal Intermittency, Phase Turbulence 
and Defect Turbulence. A detailed description can 
be found in [Chate, 1995]. 

We modify the CGLE by changing a parameter 
of the system dynamically and proportionally to the 
deviation of the system from the state to be stabi¬ 
lized. We will show that the stabilization of plane 
waves can be achieved by replacing the coefficient 
ci by ci + ^(\A\ 2 /\A pw \ 2 - 1) where 7 = 7 ,. + 17 * is 
a complex constant and \A pw \ = Aq is the modu¬ 
lus of the plane wave to be stabilized. Notice that 
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Fig. 1. Regions of the parameter space ci — C 2 for which the d = 1 CGLE displaying different kinds of regular and chaotic 
behavior. The Benjamin-Feir-Newell line (BFN line) is also shown. 


as the added term ^{\A^/Aq — 1 ) vanishes iden¬ 
tically for A — Aq, any plane wave Aq that is a 
solution of ( 1 ) is also a solution of the modified 
equation. We are not changing the solution but 
we will change its stability. The added term also 
preserves the phase invariance of the solution of the 
original CGLE, A -» Ae 1 ^, with ip being an arbi¬ 
trary phase. The modified CGLE is then explicitly 
given by, 

dtA = A + [1 + ic\ + ^{\A\ 2 /A^ — \)\d/f.A 

-{l + ic 2 )\A\ 2 A. (4) 

3. Stability Analysis 


tions for the perturbations. 

d t r = — 2 (y r k 2 + Al)r - 2Aokd x (f> — 2cikd x r 

- ciA Q d 2 x (t> + d 2 x r, (6) 

k 2 k 

d t <f> = -2c 2 A 0 r - 2yj—r - 2 C\kd x (p + 2—d x r 
Ao Aq 

T d 2 cj) + -^-d 2 r. (7) 

We consider solutions of ( 6 ) and (7) proportional to 
e pt+iqx, w j iere f or periodic boundary conditions q is 
real whereas 77 is in general a complex quantity. By 
substituting in ( 6 ) and (7) we obtain the dispersion 
relation 


We use a standard linearization procedure for 
studying the stability of the plane wave solutions 
(2) in Eq. (4). Consider the time evolution of small 
perturbations in the amplitude and phase, 


A(x, t ) = (Ao + er(x, *)) e i(**-«t+e*(z,t)) , (5) 


r] + 2(A§ + 2yrfc 2 ) + q 2 + 2ic\qk 2 iqk - c\q 2 

c\q 2 + 2 c 2 Aq - 2iqk + 2'jik 2 V + q 2 + 2ic\qk 


= 0 . 


( 8 ) 


The solutions of ( 8 ) are 


77 = —(Aq + 7 rk 2 + q 2 + 2icxqk) ±Vu + iv, (9) 


where r(x, t ) and < p(x , t) are real perturbations 
in the amplitude and phase, respectively and £ 
is a formal parameter to keep track of small 
numbers. 

Substituting (5) in (4) yields to a polynomial in 
£. The first order terms yield the linearized equa- 


where u and v are polynomials 

u — {Aq + 7 r k 2 ) 2 + 4q 2 k 2 - 2cic 2 Alq 2 - c\q A 

- 2~tiCiq 2 k 2 , ( 10 ) 

v = 4qk(ciq 2 + c 2 Aq + 7 ik 2 ). (11) 
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The real part of q indicates the growth rate of the and 

perturbations, a = (A% + 7 r k 2 ). (15) 


Re( 77 ) = -Aq - n < r k 2 



+ \fu 2 + v 2 


• ( 12 ) 


The plus and minus sign (±) correspond to two 
different branches. For the negative sign in (12), 
Re( 7 ) < 0, so these perturbations are always 
damped. 

For the positive sign of the square root in (12) 
Re(r/) = 0 at q = 0, so all the plane wave so¬ 
lutions are marginally stable. The origin of this 
neutral stability is the phase invariance A —> Ae 1 ^ 
of the solutions of Eqs. (1) and (4). For q very 
large, Re(? 7 ) ~ — q 2 , so short wavelength perturba¬ 
tions are always damped. However long wavelength 
perturbations can grow destabilizing the original 
plane wave solution, to see this we expand ( 12 ) for 
small q. 

Re(r?) = Dq 2 + 0{q A ), (13) 

where 


D= -1 


c\c 2 Al 


a 


+ 2 1 + 


C2j4q\ 


a 


+ 7i 




If this coefficient is positive, there is a range of long 
wavelength perturbations that grow. The condition 
D < 0 is necessary for stability but not sufficient, 
since the growth coefficient obtained from the full 
expression ( 12 ) can be positive for some finite value 
of q despite the coefficient D being negative. How¬ 
ever, for the values of c\ and c 2 considered in this 
work the requirement D < 0 gives a very good limit 
for the stability region. 

For the unperturbed CGLE (7 = 0) the condi¬ 
tion D < 0, leads to the standard Eckhaus instabil¬ 
ity limit: |fc| < kE with kE given by (3). Since the 
control is done through a diffusive terhi, the added 
term never changes the stability of the homogeneous 
solution k = 0. For an arbitrary plane wave of wave 
number k one has to solve the cubic equation D = 0 
to find explicitly the range of values of k for which 
the plane wave is stable. In Figs. 2 and 3 we plot 
this range as a function of the parameter c\ for sev¬ 
eral values of 7 and c 2 . When no control is applied 
(7 = 0) the stability region in the k-c\ plane is lim¬ 
ited by a branch of Eq. (3) (dashed line in Figs. 2 
and 3) whose vertex corresponds to the BFN point. 
Decreasing the value of ci the width of the stabil¬ 
ity region |A;| < kE increases, and for c\ -> — 00, 
kE -+ 1 . 



Ikl Ikl 


Fig. 2. The shadowed zone shows the stability region for the 
plane wave (2) for ci = —0.9. Left column and from top to 
bottom, 7 i = 0 and = 1, 2, 3. Right column and from 
top to bottom, 7 i = 1 and 7 = 1, 2, 3. For comparison, the 
boundary of the stability region for 7 = 0, given by Eq. (3), 
is shown in all the figures as a dashed line. 
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Fig. 3. The shadowed zone shows the stability region for the 
plane wave (2) for C 2 = —2.1. Left column and from top to 
bottom, 7 i = 0 and y r = 1, 2, 3. Right column and from 
top to bottom, ji = 1 and 7 r = 1, 2, 3. For comparison, the 
boundary of the stability region for 7 = 0 is shown in all the 
figures as a dashed line. 
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As |'y| is increased the stability region changes 
its shape, so that it is possible to stabilize plane 
waves which were unstable without control. In¬ 
creasing the value of |y| one finds stabilization of 
plane waves for values of c\ well above the BFN 
line, in the regions of phase turbulence and defect 
turbulence of the original CGLE. In particular, for 
C 2 = —0.9 the BFN line corresponds in Fig. 2 to 
a horizontal line at ci = 1.11. Below this value 
for ci, plane waves with a large wavevector (out¬ 


side the dashed line) are unstable without control. 
As shown in Fig. 2 they can be easily stabilized 
with a real control parameter. For ci > 1.11 we 
are in the phase turbulence regime up to ci « 2.3 
(see Fig. 1). For larger values of ci we enter in the 
defect turbulence region. It is possible to stabilize 
plane waves in the region of phase turbulence tak¬ 
ing 7 ,. > 3 and 71 = 0, or 7 r > 1 and 7 * = 1. For 
C 2 = —2.1 the BFN line corresponds in Fig. 3 to 
a horizontal line at ci = 0.476. Below that level 




Fig. 4. Spatiotemporal evolution of the CGLE (4) for ci = 1.5, C 2 = —0.9 starting from a perturbed plane wave (16) with 
k = 0.5 and cr — 0.007. Figures (a) and (b) show \A(x, t)| with time running upwards from t = 0 to t = 1000 and x in the 
horizontal direction for 7 = 0 and 7 = 3 + Oi, respectively. The absolute value of the field \A(x, to)| and the phase gradient 
d x 4>(x, to) at to = 950 are displayed in (c) and (d) for 7 = 0 and 7 = 3 + Oi, respectively. 








1854 R. Montague <h P. Colet 


there are the coexistence of stable plane waves 
and spatiotemporal intermitency. Plane waves with 
wavevector outside the dashed line are unstable and 
will lead to spatiotemporal intermitency without 
control. As shown in Fig. 3 the real part of 7 is very 
effective in the stabilization of these plane waves. 
For c\ > 0.476 we are above the BFN line, and in¬ 
creasing ci we enter first in the bichaos regime and 
later in the defect turbulence regime (see Fig. 1). 
As shown in Fig. 3, stabilization of plane waves in 
both regions is possible using for example 7 = 1 + i. 


In general, the real part of 7 is specially effec¬ 
tive in extending the stability region to plane waves 
of a large wave number. We should stress that sta¬ 
bilization of plane waves with wavevector close to 
one was extremely difficult to achieve using only 
an imaginary nonlinear diffusion term, even tak¬ 
ing large values for 7 j, as shown in [Montagne & 
Colet, 1997]. Above the BFN line, stabilization of 
unstable plane waves in the regions of phase and 
defect turbulence can be achieved taking 7 to be 
purely real (left column in Figs. 2 and 3) or purely 
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Fig. 5. Spatiotemporal evolution of the CGLE for ci = 2.5, C 2 = —2.1 starting from a perturbed plane wave with k = 0.9. 
Figures (a) and (b) show |4(a:, t)| for 7 = 0 and 7 = 1 + i, respectively. The values of \A(x, to)| and d x (j>{x, t) at to = 950 are 
displayed in (c) and (d) for 7 = 0 and 7=1 + 2 , respectively. Other parameters are as in Fig. 4. 
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imaginary (see ref. [Montagne & Colet, 1997]). 
However, it is more efficient to take 7 complex, since 
one can achieve stabilization with smaller values 
of | 7 |. 

4. Numerical simulations 

We have integrated numerically Eqs. (1) and (4) us¬ 
ing a pseudospectral code with periodic boundary 
conditions and second-order accuracy in time. Spa¬ 
tial resolution was typically 1024 modes. Time step 
was typically At = 0.001. The system size was al¬ 
ways taken as L = 512. The details of the numerical 
method can be seen in [Montagne et al., 1997]. We 
start from an initial condition which corresponds to 
a plane wave plus a small random perturbation 

A(x, t = 0 ) = Vl — k 2 e lkx + &i{x) (16) 

where £(x) is a complex Gaussian random pertur¬ 
bation of zero mean and variance (£(x)£* (x 1 )) — 
2S(x — x'). 

We have performed numerical simulations in 
different regions of the phase diagram (Fig. 1) to 
verify the results obtained from the linear stability 
analysis when finite size perturbations are applied. 
We have found a very good agreement between the 
prediction of the linear stability analysis and the 
numerical simulations. As characteristic examples 
we show the following results. 

Figure 4 shows the stabilization of a plane wave 
for parameter values c\ = 1.5 and C 2 = —0.9. This 
corresponds to the phase turbulence regime (see 
Fig. 1), where no plane waves are stable for the 
original CGLE. As predicted by the linear stability 
analysis, a perturbed plane wave with k = 0.5 can 
easily be stabilized with j r = 3 and 7 , = 0 while 
for 7 = 0 the same initial condition decays in time 
t = 80 (approx.) to phase turbulence. 

Figure 5, obtained for parameter values c\ — 
2.5 and C 2 = —2.1, shows stabilization of plane 
waves in the region of defect turbulence, where for 
the unperturbed CGLE there are no stable plane 
waves and the field A shows a strongly disordered 
STC state characterized by the presence of defects. 
As predicted by the linear stability analysis, a per¬ 
turbed plane wave with k = 0.9 can be stabilized 
with 7 r = 1 and 7 ■{ = 1 . 

5. Concluding Remarks 

Our study of the dynamics of the CGLE with a 
nonlinear diffusion term demonstrates the possibil¬ 


ity of stabilization of an unstable plane wave in dif¬ 
ferent regimes of STC. Since the added term van¬ 
ish for the stabilized plane wave, this plane wave 
is exactly the same unstable solution of the origi¬ 
nal CGLE. Although our method does not change 
the stability of the homogeneous solution k = 0, 
it is quite effective in stabilizing plane waves with 
nonzero wavevector. We have shown that the real 
part of the nonlinear diffusive term is specially ef¬ 
fective for stabilizing plane waves with wavevector 
close to one. Also, in general it is more efficient to 
use a nonlinear diffusive term with both real and 
imaginary parts, in the sense that stabilization can 
be achieved with a control parameter with smaller 
absolute value. 

Finally, we have studied numerically the stabil¬ 
ity of the plane waves when finite size perturbations 
are applied. The results are in excellent agreement 
with the analytical predictions of the linear stability 
analysis. 
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